Context Navigation

spacecurve_mod.F90 @ 141

Revision 141, 50.6 KB checked in by maltrud, 4 years ago (diff)

LANL merge with NCAR revision 17212. LANL version removed several features,
including overflows, run-time bathymetry generation, and some of the
coupling infrastructure. Also, some of the tavg calls are different, which
will be addressed in a later update. High frequency pointwise and regional
diagnostics have also been removed, to be replaced in the future by the "moorings"
infrastructure. Diagnostics using "lat_aux_grid" have also been removed.

Line
1	!BOP
2	! !MODULE: spacecurve_mod
3
4	module spacecurve_mod
5
6	! !DESCRIPTION:
7	! This module contains routines necessary to
8	! create space-filling curves.
9	!
10	! !REVISION HISTORY:
11	! CVS: $Id: spacecurve_mod.F90,v 1.3 2006/07/06 15:48:09 dennis Exp $
12	! CVS: $Name: $
13
14	! !USES:
15	use kinds_mod
16
17	implicit none
18
19	! !PUBLIC TYPES:
20
21	type, public :: factor_t
22	integer(int_kind) :: numfact ! The # of factors for a value
23	integer(int_kind), dimension(:), pointer :: factors ! The factors
24	integer(int_kind), dimension(:), pointer :: used
25	end type
26
27	! !PUBLIC MEMBER FUNCTIONS:
28
29	public :: GenSpaceCurve, &
30	IsLoadBalanced
31
32	public :: Factor, &
33	IsFactorable, &
34	PrintFactor, &
35	ProdFactor, &
36	MatchFactor
37
38	! !PRIVATE MEMBER FUNCTIONS:
39
40	private :: map, &
41	PeanoM, &
42	Hilbert, &
43	Cinco, &
44	GenCurve
45
46	private :: FirstFactor, &
47	FindandMark
48
49	integer(int_kind), dimension(:,:), allocatable :: &
50	dir, &! direction to move along each level
51	ordered, &! the ordering
52	traversal ! index increments to traverse domain
53
54	integer(int_kind), dimension(:), allocatable :: &
55	pos ! position along each of the axes
56
57	integer(int_kind) :: &
58	maxdim, &! dimensionality of entire space
59	vcnt ! visitation count
60
61	logical :: verbose=.FALSE.
62
63	type (factor_t), public,save :: fact ! stores the factorization
64
65	!EOP
66	!EOC
67	!***********************************************************************
68
69	contains
70
71	!***********************************************************************
72	!BOP
73	! !IROUTINE: Cinco
74	! !INTERFACE:
75
76	recursive function Cinco(l,type,ma,md,ja,jd) result(ierr)
77
78	! !DESCRIPTION:
79	! This subroutine implements a Cinco space-filling curve.
80	! Cinco curves connect a Nb x Nb block of points where
81	!
82	! Nb = 5^p
83	!
84	! !REVISION HISTORY:
85	! same as module
86	!
87
88
89	! !INPUT PARAMETERS
90
91	integer(int_kind), intent(in) :: &
92	l, & ! level of the space-filling curve
93	type, & ! type of SFC curve
94	ma, & ! Major axis [0,1]
95	md, & ! direction of major axis [-1,1]
96	ja, & ! joiner axis [0,1]
97	jd ! direction of joiner axis [-1,1]
98
99	! !OUTPUT PARAMETERS
100
101	integer(int_kind) :: ierr ! error return code
102
103	!EOP
104	!BOC
105	!-----------------------------------------------------------------------
106	!
107	! local variables
108	!
109	!-----------------------------------------------------------------------
110
111	integer(int_kind) :: &
112	lma, &! local major axis (next level)
113	lmd, &! local major direction (next level)
114	lja, &! local joiner axis (next level)
115	ljd, &! local joiner direction (next level)
116	ltype, &! type of SFC on next level
117	ll ! next level down
118
119	logical :: debug = .FALSE.
120
121	!-----------------------------------------------------------------------
122	ll = l
123	if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
124
125	!--------------------------------------------------------------
126	! Position [0,0]
127	!--------------------------------------------------------------
128	lma = ma
129	lmd = md
130	lja = lma
131	ljd = lmd
132
133	if(ll .gt. 1) then
134	if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
135	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
136	if(debug) call PrintCurve(ordered)
137	else
138	ierr = IncrementCurve(lja,ljd)
139	if(debug) write(,) 'Cinco: After Position [0,0] ',pos
140	endif
141
142	!--------------------------------------------------------------
143	! Position [1,0]
144	!--------------------------------------------------------------
145	lma = ma
146	lmd = md
147	lja = lma
148	ljd = lmd
149
150	if(ll .gt. 1) then
151	if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
152	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
153	if(debug) call PrintCurve(ordered)
154	else
155	ierr = IncrementCurve(lja,ljd)
156	if(debug) write(,) 'After Position [1,0] ',pos
157	endif
158
159	!--------------------------------------------------------------
160	! Position [2,0]
161	!--------------------------------------------------------------
162	lma = MOD(ma+1,maxdim)
163	lmd = md
164	lja = lma
165	ljd = lmd
166
167	if(ll .gt. 1) then
168	if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
169	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
170	if(debug) call PrintCurve(ordered)
171	else
172	ierr = IncrementCurve(lja,ljd)
173	if(debug) write(,) 'After Position [2,0] ',pos
174	endif
175
176	!--------------------------------------------------------------
177	! Position [2,1]
178	!--------------------------------------------------------------
179	lma = MOD(ma+1,maxdim)
180	lmd = md
181	lja = lma
182	ljd = lmd
183
184	if(ll .gt. 1) then
185	if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
186	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
187	if(debug) call PrintCurve(ordered)
188	else
189	ierr = IncrementCurve(lja,ljd)
190	if(debug) write(,) 'After Position [2,1] ',pos
191	endif
192
193	!--------------------------------------------------------------
194	! Position [2,2]
195	!--------------------------------------------------------------
196	lma = MOD(ma+1,maxdim)
197	lmd = md
198	lja = ma
199	ljd = -md
200
201	if(ll .gt. 1) then
202	if(debug) write(*,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
203	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
204	if(debug) call PrintCurve(ordered)
205	else
206	ierr = IncrementCurve(lja,ljd)
207	if(debug) write(,) 'After Position [2,2] ',pos
208	endif
209
210	!--------------------------------------------------------------
211	! Position [1,2]
212	!--------------------------------------------------------------
213	lma = MOD(ma+1,maxdim)
214	lmd = -md
215	lja = lma
216	ljd = lmd
217
218	if(ll .gt. 1) then
219	if(debug) write(*,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
220	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
221	if(debug) call PrintCurve(ordered)
222	else
223	ierr = IncrementCurve(lja,ljd)
224	if(debug) write(,) 'After Position [1,2] ',pos
225	endif
226
227	!--------------------------------------------------------------
228	! Position [1,1]
229	!--------------------------------------------------------------
230	lma = ma
231	lmd = -md
232	lja = lma
233	ljd = lmd
234
235	if(ll .gt. 1) then
236	if(debug) write(*,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
237	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
238	if(debug) call PrintCurve(ordered)
239	else
240	ierr = IncrementCurve(lja,ljd)
241	if(debug) write(,) 'After Position [1,1] ',pos
242	endif
243
244	!--------------------------------------------------------------
245	! Position [0,1]
246	!--------------------------------------------------------------
247	lma = ma
248	lmd = -md
249	lja = MOD(ma+1,maxdim)
250	ljd = md
251
252	if(ll .gt. 1) then
253	if(debug) write(*,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
254	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
255	if(debug) call PrintCurve(ordered)
256	else
257	ierr = IncrementCurve(lja,ljd)
258	if(debug) write(,) 'After Position [0,1] ',pos
259	endif
260
261	!--------------------------------------------------------------
262	! Position [0,2]
263	!--------------------------------------------------------------
264	lma = MOD(ma+1,maxdim)
265	lmd = md
266	lja = lma
267	ljd = lmd
268
269	if(ll .gt. 1) then
270	if(debug) write(*,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
271	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
272	if(debug) call PrintCurve(ordered)
273	else
274	ierr = IncrementCurve(lja,ljd)
275	if(debug) write(,) 'After Position [0,2] ',pos
276	endif
277
278	!--------------------------------------------------------------
279	! Position [0,3]
280	!--------------------------------------------------------------
281	lma = MOD(ma+1,maxdim)
282	lmd = md
283	lja = lma
284	ljd = lmd
285
286	if(ll .gt. 1) then
287	if(debug) write(*,30) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
288	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
289	if(debug) call PrintCurve(ordered)
290	else
291	ierr = IncrementCurve(lja,ljd)
292	if(debug) write(,) 'After Position [0,3] ',pos
293	endif
294
295	!--------------------------------------------------------------
296	! Position [0,4]
297	!--------------------------------------------------------------
298	lma = ma
299	lmd = md
300	lja = lma
301	ljd = lmd
302
303	if(ll .gt. 1) then
304	if(debug) write(*,31) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
305	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
306	if(debug) call PrintCurve(ordered)
307	else
308	ierr = IncrementCurve(lja,ljd)
309	if(debug) write(,) 'After Position [0,4] ',pos
310	endif
311
312	!--------------------------------------------------------------
313	! Position [1,4]
314	!--------------------------------------------------------------
315	lma = ma
316	lmd = md
317	lja = MOD(ma+1,maxdim)
318	ljd = -md
319
320	if(ll .gt. 1) then
321	if(debug) write(*,32) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
322	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
323	if(debug) call PrintCurve(ordered)
324	else
325	ierr = IncrementCurve(lja,ljd)
326	if(debug) write(,) 'After Position [1,4] ',pos
327	endif
328
329	!--------------------------------------------------------------
330	! Position [1,3]
331	!--------------------------------------------------------------
332	lma = MOD(ma+1,maxdim)
333	lmd = -md
334	lja = ma
335	ljd = md
336
337	if(ll .gt. 1) then
338	if(debug) write(*,33) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
339	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
340	if(debug) call PrintCurve(ordered)
341	else
342	ierr = IncrementCurve(lja,ljd)
343	if(debug) write(,) 'After Position [1,3] ',pos
344	endif
345
346	!--------------------------------------------------------------
347	! Position [2,3]
348	!--------------------------------------------------------------
349	lma = MOD(ma+1,maxdim)
350	lmd = md
351	lja = lma
352	ljd = lmd
353
354	if(ll .gt. 1) then
355	if(debug) write(*,34) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
356	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
357	if(debug) call PrintCurve(ordered)
358	else
359	ierr = IncrementCurve(lja,ljd)
360	if(debug) write(,) 'After Position [2,3] ',pos
361	endif
362
363	!--------------------------------------------------------------
364	! Position [2,4]
365	!--------------------------------------------------------------
366	lma = ma
367	lmd = md
368	lja = lma
369	ljd = lmd
370
371	if(ll .gt. 1) then
372	if(debug) write(*,35) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
373	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
374	if(debug) call PrintCurve(ordered)
375	else
376	ierr = IncrementCurve(lja,ljd)
377	if(debug) write(,) 'After Position [2,4] ',pos
378	endif
379
380	!--------------------------------------------------------------
381	! Position [3,4]
382	!--------------------------------------------------------------
383	lma = ma
384	lmd = md
385	lja = lma
386	ljd = lmd
387
388	if(ll .gt. 1) then
389	if(debug) write(*,36) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
390	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
391	if(debug) call PrintCurve(ordered)
392	else
393	ierr = IncrementCurve(lja,ljd)
394	if(debug) write(,) 'After Position [3,4] ',pos
395	endif
396
397	!--------------------------------------------------------------
398	! Position [4,4]
399	!--------------------------------------------------------------
400	lma = ma
401	lmd = md
402	lja = MOD(ma+1,maxdim)
403	ljd = -md
404
405	if(ll .gt. 1) then
406	if(debug) write(*,37) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
407	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
408	if(debug) call PrintCurve(ordered)
409	else
410	ierr = IncrementCurve(lja,ljd)
411	if(debug) write(,) 'After Position [4,4] ',pos
412	endif
413
414	!--------------------------------------------------------------
415	! Position [4,3]
416	!--------------------------------------------------------------
417	lma = ma
418	lmd = -md
419	lja = lma
420	ljd = lmd
421
422	if(ll .gt. 1) then
423	if(debug) write(*,38) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
424	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
425	if(debug) call PrintCurve(ordered)
426	else
427	ierr = IncrementCurve(lja,ljd)
428	if(debug) write(,) 'After Position [4,3] ',pos
429	endif
430
431	!--------------------------------------------------------------
432	! Position [3,3]
433	!--------------------------------------------------------------
434	lma = MOD(ma+1,maxdim)
435	lmd = -md
436	lja = lma
437	ljd = lmd
438
439	if(ll .gt. 1) then
440	if(debug) write(*,39) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
441	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
442	if(debug) call PrintCurve(ordered)
443	else
444	ierr = IncrementCurve(lja,ljd)
445	if(debug) write(,) 'After Position [3,3] ',pos
446	endif
447
448	!--------------------------------------------------------------
449	! Position [3,2]
450	!--------------------------------------------------------------
451	lma = MOD(ma+1,maxdim)
452	lmd = -md
453	lja = ma
454	ljd = md
455
456	if(ll .gt. 1) then
457	if(debug) write(*,40) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
458	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
459	if(debug) call PrintCurve(ordered)
460	else
461	ierr = IncrementCurve(lja,ljd)
462	if(debug) write(,) 'After Position [3,2] ',pos
463	endif
464
465	!--------------------------------------------------------------
466	! Position [4,2]
467	!--------------------------------------------------------------
468	lma = ma
469	lmd = md
470	lja = MOD(ma+1,maxdim)
471	ljd = -md
472
473	if(ll .gt. 1) then
474	if(debug) write(*,41) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
475	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
476	if(debug) call PrintCurve(ordered)
477	else
478	ierr = IncrementCurve(lja,ljd)
479	if(debug) write(,) 'After Position [4,2] ',pos
480	endif
481
482	!--------------------------------------------------------------
483	! Position [4,1]
484	!--------------------------------------------------------------
485	lma = ma
486	lmd = -md
487	lja = lma
488	ljd = lmd
489
490	if(ll .gt. 1) then
491	if(debug) write(*,42) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
492	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
493	if(debug) call PrintCurve(ordered)
494	else
495	ierr = IncrementCurve(lja,ljd)
496	if(debug) write(,) 'After Position [4,1] ',pos
497	endif
498
499	!--------------------------------------------------------------
500	! Position [3,1]
501	!--------------------------------------------------------------
502	lma = MOD(ma+1,maxdim)
503	lmd = -md
504	lja = lma
505	ljd = lmd
506
507	if(ll .gt. 1) then
508	if(debug) write(*,43) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
509	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
510	if(debug) call PrintCurve(ordered)
511	else
512	ierr = IncrementCurve(lja,ljd)
513	if(debug) write(,) 'After Position [3,1] ',pos
514	endif
515
516	!--------------------------------------------------------------
517	! Position [3,0]
518	!--------------------------------------------------------------
519	lma = MOD(ma+1,maxdim)
520	lmd = -md
521	lja = ma
522	ljd = md
523
524	if(ll .gt. 1) then
525	if(debug) write(*,44) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
526	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
527	if(debug) call PrintCurve(ordered)
528	else
529	ierr = IncrementCurve(lja,ljd)
530	if(debug) write(,) 'After Position [3,0] ',pos
531	endif
532
533	!--------------------------------------------------------------
534	! Position [4,0]
535	!--------------------------------------------------------------
536	lma = ma
537	lmd = md
538	lja = ja
539	ljd = jd
540
541	if(ll .gt. 1) then
542	if(debug) write(*,45) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
543	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
544	if(debug) call PrintCurve(ordered)
545	else
546	ierr = IncrementCurve(lja,ljd)
547	if(debug) write(,) 'After Position [4,0] ',pos
548	endif
549
550	21 format('Call Cinco Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
551	22 format('Call Cinco Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
552	23 format('Call Cinco Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
553	24 format('Call Cinco Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
554	25 format('Call Cinco Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
555	26 format('Call Cinco Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
556	27 format('Call Cinco Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
557	28 format('Call Cinco Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
558	29 format('Call Cinco Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
559	30 format('Call Cinco Pos [0,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
560	31 format('Call Cinco Pos [0,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
561	32 format('Call Cinco Pos [1,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
562	33 format('Call Cinco Pos [1,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
563	34 format('Call Cinco Pos [2,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
564	35 format('Call Cinco Pos [2,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
565	36 format('Call Cinco Pos [3,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
566	37 format('Call Cinco Pos [4,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
567	38 format('Call Cinco Pos [4,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
568	39 format('Call Cinco Pos [3,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
569	40 format('Call Cinco Pos [3,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
570	41 format('Call Cinco Pos [4,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
571	42 format('Call Cinco Pos [4,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
572	43 format('Call Cinco Pos [3,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
573	44 format('Call Cinco Pos [3,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
574	45 format('Call Cinco Pos [4,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
575
576	!EOC
577	!-----------------------------------------------------------------------
578
579	end function Cinco
580
581	!***********************************************************************
582	!BOP
583	! !IROUTINE: PeanoM
584	! !INTERFACE:
585
586	recursive function PeanoM(l,type,ma,md,ja,jd) result(ierr)
587
588	! !DESCRIPTION:
589	! This function implements a meandering Peano
590	! space-filling curve. A meandering Peano curve
591	! connects a Nb x Nb block of points where
592	!
593	! Nb = 3^p
594	!
595	! !REVISION HISTORY:
596	! same as module
597	!
598
599	! !INPUT PARAMETERS
600
601	integer(int_kind), intent(in) :: &
602	l, & ! level of the space-filling curve
603	type, & ! type of SFC curve
604	ma, & ! Major axis [0,1]
605	md, & ! direction of major axis [-1,1]
606	ja, & ! joiner axis [0,1]
607	jd ! direction of joiner axis [-1,1]
608
609	! !OUTPUT PARAMETERS
610
611	integer(int_kind) :: ierr ! error return code
612
613	!EOP
614	!BOC
615	!-----------------------------------------------------------------------
616	!
617	! local variables
618	!
619	!-----------------------------------------------------------------------
620
621
622	integer(int_kind) :: &
623	lma, &! local major axis (next level)
624	lmd, &! local major direction (next level)
625	lja, &! local joiner axis (next level)
626	ljd, &! local joiner direction (next level)
627	ltype, &! type of SFC on next level
628	ll ! next level down
629
630	logical :: debug = .FALSE.
631
632	!-----------------------------------------------------------------------
633
634	ll = l
635	if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
636	!--------------------------------------------------------------
637	! Position [0,0]
638	!--------------------------------------------------------------
639	lma = MOD(ma+1,maxdim)
640	lmd = md
641	lja = lma
642	ljd = lmd
643
644	if(ll .gt. 1) then
645	if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
646	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
647	if(debug) call PrintCurve(ordered)
648	else
649	ierr = IncrementCurve(lja,ljd)
650	if(debug) write(,) 'PeanoM: After Position [0,0] ',pos
651	endif
652
653
654	!--------------------------------------------------------------
655	! Position [0,1]
656	!--------------------------------------------------------------
657	lma = MOD(ma+1,maxdim)
658	lmd = md
659	lja = lma
660	ljd = lmd
661	if(ll .gt. 1) then
662	if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
663	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
664	if(debug) call PrintCurve(ordered)
665	else
666	ierr = IncrementCurve(lja,ljd)
667	if(debug) write(,) 'PeanoM: After Position [0,1] ',pos
668	endif
669
670	!--------------------------------------------------------------
671	! Position [0,2]
672	!--------------------------------------------------------------
673	lma = ma
674	lmd = md
675	lja = lma
676	ljd = lmd
677	if(ll .gt. 1) then
678	if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
679	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
680	if(debug) call PrintCurve(ordered)
681	else
682	ierr = IncrementCurve(lja,ljd)
683	if(debug) write(,) 'PeanoM: After Position [0,2] ',pos
684	endif
685
686	!--------------------------------------------------------------
687	! Position [1,2]
688	!--------------------------------------------------------------
689	lma = ma
690	lmd = md
691	lja = lma
692	ljd = lmd
693	if(ll .gt. 1) then
694	if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
695	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
696	if(debug) call PrintCurve(ordered)
697	else
698	ierr = IncrementCurve(lja,ljd)
699	if(debug) write(,) 'PeanoM: After Position [1,2] ',pos
700	endif
701
702
703	!--------------------------------------------------------------
704	! Position [2,2]
705	!--------------------------------------------------------------
706	lma = ma
707	lmd = md
708	lja = MOD(lma+1,maxdim)
709	ljd = -lmd
710
711	if(ll .gt. 1) then
712	if(debug) write(*,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
713	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
714	if(debug) call PrintCurve(ordered)
715	else
716	ierr = IncrementCurve(lja,ljd)
717	if(debug) write(,) 'PeanoM: After Position [2,2] ',pos
718	endif
719
720	!--------------------------------------------------------------
721	! Position [2,1]
722	!--------------------------------------------------------------
723	lma = ma
724	lmd = -md
725	lja = lma
726	ljd = lmd
727
728	if(ll .gt. 1) then
729	if(debug) write(*,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
730	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
731	if(debug) call PrintCurve(ordered)
732	else
733	ierr = IncrementCurve(lja,ljd)
734	if(debug) write(,) 'PeanoM: After Position [2,1] ',pos
735	endif
736
737	!--------------------------------------------------------------
738	! Position [1,1]
739	!--------------------------------------------------------------
740	lma = MOD(ma+1,maxdim)
741	lmd = -md
742	lja = lma
743	ljd = lmd
744
745	if(ll .gt. 1) then
746	if(debug) write(*,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
747	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
748	if(debug) call PrintCurve(ordered)
749	else
750	ierr = IncrementCurve(lja,ljd)
751	if(debug) write(,) 'PeanoM: After Position [1,1] ',pos
752	endif
753
754
755	!--------------------------------------------------------------
756	! Position [1,0]
757	!--------------------------------------------------------------
758	lma = MOD(ma+1,maxdim)
759	lmd = -md
760	lja = MOD(lma+1,maxdim)
761	ljd = -lmd
762
763	if(ll .gt. 1) then
764	if(debug) write(*,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
765	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
766	if(debug) call PrintCurve(ordered)
767	else
768	ierr = IncrementCurve(lja,ljd)
769	if(debug) write(,) 'PeanoM: After Position [1,0] ',pos
770	endif
771
772	!--------------------------------------------------------------
773	! Position [2,0]
774	!--------------------------------------------------------------
775	lma = ma
776	lmd = md
777	lja = ja
778	ljd = jd
779
780	if(ll .gt. 1) then
781	if(debug) write(*,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
782	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
783	if(debug) call PrintCurve(ordered)
784	else
785	ierr = IncrementCurve(lja,ljd)
786	if(debug) write(,) 'PeanoM: After Position [2,0] ',pos
787	endif
788
789	21 format('Call PeanoM Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
790	22 format('Call PeanoM Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
791	23 format('Call PeanoM Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
792	24 format('Call PeanoM Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
793	25 format('Call PeanoM Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
794	26 format('Call PeanoM Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
795	27 format('Call PeanoM Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
796	28 format('Call PeanoM Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
797	29 format('Call PeanoM Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
798
799	!EOC
800	!-----------------------------------------------------------------------
801
802	end function PeanoM
803
804	!***********************************************************************
805	!BOP
806	! !IROUTINE: Hilbert
807	! !INTERFACE:
808
809	recursive function Hilbert(l,type,ma,md,ja,jd) result(ierr)
810
811	! !DESCRIPTION:
812	! This function implements a Hilbert space-filling curve.
813	! A Hilbert curve connect a Nb x Nb block of points where
814	!
815	! Nb = 2^p
816	!
817	! !REVISION HISTORY:
818	! same as module
819	!
820
821
822	! !INPUT PARAMETERS
823
824	integer(int_kind), intent(in) :: &
825	l, & ! level of the space-filling curve
826	type, & ! type of SFC curve
827	ma, & ! Major axis [0,1]
828	md, & ! direction of major axis [-1,1]
829	ja, & ! joiner axis [0,1]
830	jd ! direction of joiner axis [-1,1]
831
832	! !OUTPUT PARAMETERS
833
834	integer(int_kind) :: ierr ! error return code
835
836	!EOP
837	!BOC
838	!-----------------------------------------------------------------------
839	!
840	! local variables
841	!
842	!-----------------------------------------------------------------------
843
844
845	integer(int_kind) :: &
846	lma, &! local major axis (next level)
847	lmd, &! local major direction (next level)
848	lja, &! local joiner axis (next level)
849	ljd, &! local joiner direction (next level)
850	ltype, &! type of SFC on next level
851	ll ! next level down
852
853	logical :: debug = .FALSE.
854
855	!-----------------------------------------------------------------------
856	ll = l
857	if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
858	!--------------------------------------------------------------
859	! Position [0,0]
860	!--------------------------------------------------------------
861	lma = MOD(ma+1,maxdim)
862	lmd = md
863	lja = lma
864	ljd = lmd
865
866	if(ll .gt. 1) then
867	if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
868	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
869	if(debug) call PrintCurve(ordered)
870	else
871	ierr = IncrementCurve(lja,ljd)
872	if(debug) write(,) 'Hilbert: After Position [0,0] ',pos
873	endif
874
875
876	!--------------------------------------------------------------
877	! Position [0,1]
878	!--------------------------------------------------------------
879	lma = ma
880	lmd = md
881	lja = lma
882	ljd = lmd
883	if(ll .gt. 1) then
884	if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
885	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
886	if(debug) call PrintCurve(ordered)
887	else
888	ierr = IncrementCurve(lja,ljd)
889	if(debug) write(,) 'Hilbert: After Position [0,1] ',pos
890	endif
891
892
893	!--------------------------------------------------------------
894	! Position [1,1]
895	!--------------------------------------------------------------
896	lma = ma
897	lmd = md
898	lja = MOD(ma+1,maxdim)
899	ljd = -md
900
901	if(ll .gt. 1) then
902	if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
903	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
904	if(debug) call PrintCurve(ordered)
905	else
906	ierr = IncrementCurve(lja,ljd)
907	if(debug) write(,) 'Hilbert: After Position [1,1] ',pos
908	endif
909
910	!--------------------------------------------------------------
911	! Position [1,0]
912	!--------------------------------------------------------------
913	lma = MOD(ma+1,maxdim)
914	lmd = -md
915	lja = ja
916	ljd = jd
917
918	if(ll .gt. 1) then
919	if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
920	ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
921	if(debug) call PrintCurve(ordered)
922	else
923	ierr = IncrementCurve(lja,ljd)
924	if(debug) write(,) 'Hilbert: After Position [1,0] ',pos
925	endif
926
927	21 format('Call Hilbert Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
928	22 format('Call Hilbert Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
929	23 format('Call Hilbert Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
930	24 format('Call Hilbert Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
931
932	!EOC
933	!-----------------------------------------------------------------------
934
935	end function hilbert
936
937	!***********************************************************************
938	!BOP
939	! !IROUTINE: IncrementCurve
940	! !INTERFACE:
941
942	function IncrementCurve(ja,jd) result(ierr)
943
944	! !DESCRIPTION:
945	! This function creates the curve which is store in the
946	! the ordered array. The curve is implemented by
947	! incrementing the curve in the direction [jd] of axis [ja].
948	!
949	! !REVISION HISTORY:
950	! same as module
951	!
952
953	! !INPUT PARAMETERS:
954	integer(int_kind) :: &
955	ja, &! axis to increment
956	jd ! direction along axis
957
958	! !OUTPUT PARAMETERS:
959	integer(int_kind) :: ierr ! error return code
960
961	!-----------------------------
962	! mark the newly visited point
963	!-----------------------------
964	ordered(pos(0)+1,pos(1)+1) = vcnt
965
966	!------------------------------------
967	! increment curve and update position
968	!------------------------------------
969	vcnt = vcnt + 1
970	pos(ja) = pos(ja) + jd
971
972	ierr = 0
973	!EOC
974	!-----------------------------------------------------------------------
975
976	end function IncrementCurve
977
978	!***********************************************************************
979	!BOP
980	! !IROUTINE: log2
981	! !INTERFACE:
982
983	function log2( n)
984
985	! !DESCRIPTION:
986	! This function calculates the log2 of its integer
987	! input.
988	!
989	! !REVISION HISTORY:
990	! same as module
991
992	! !INPUT PARAMETERS:
993
994	integer(int_kind), intent(in) :: n ! integer value to find the log2
995
996	! !OUTPUT PARAMETERS:
997
998	integer(int_kind) :: log2
999
1000	!EOP
1001	!BOC
1002	!-----------------------------------------------------------------------
1003	!
1004	! local variables
1005	!
1006	!-----------------------------------------------------------------------
1007
1008	integer(int_kind) :: tmp
1009
1010	!-------------------------------
1011	! Find the log2 of input value
1012	!-------------------------------
1013	log2 = 1
1014	tmp =n
1015	do while (tmp/2 .ne. 1)
1016	tmp=tmp/2
1017	log2=log2+1
1018	enddo
1019
1020	!EOP
1021	!-----------------------------------------------------------------------
1022
1023	end function log2
1024
1025	!***********************************************************************
1026	!BOP
1027	! !IROUTINE: IsLoadBalanced
1028	! !INTERFACE:
1029
1030	function IsLoadBalanced(nelem,npart)
1031
1032	! !DESCRIPTION:
1033	! This function determines if we can create
1034	! a perfectly load-balanced partitioning.
1035	!
1036	! !REVISION HISTORY:
1037	! same as module
1038
1039	! !INTPUT PARAMETERS:
1040
1041	integer(int_kind), intent(in) :: &
1042	nelem, & ! number of blocks/elements to partition
1043	npart ! size of partition
1044
1045	! !OUTPUT PARAMETERS:
1046	logical :: IsLoadBalanced ! .TRUE. if a perfectly load balanced
1047	! partition is possible
1048	!EOP
1049	!BOC
1050	!-----------------------------------------------------------------------
1051	!
1052	! local variables
1053	!
1054	!-----------------------------------------------------------------------
1055
1056	integer(int_kind) :: tmp1 ! temporary int
1057
1058	!-----------------------------------------------------------------------
1059	tmp1 = nelem/npart
1060
1061	if(npart*tmp1 == nelem ) then
1062	IsLoadBalanced=.TRUE.
1063	else
1064	IsLoadBalanced=.FALSE.
1065	endif
1066
1067	!EOP
1068	!-----------------------------------------------------------------------
1069
1070	end function IsLoadBalanced
1071
1072	!***********************************************************************
1073	!BOP
1074	! !IROUTINE: GenCurve
1075	! !INTERFACE:
1076
1077	function GenCurve(l,type,ma,md,ja,jd) result(ierr)
1078
1079	! !DESCRIPTION:
1080	! This subroutine generates the next level down
1081	! space-filling curve
1082	!
1083	! !REVISION HISTORY:
1084	! same as module
1085	!
1086
1087	! !INPUT PARAMETERS
1088
1089	integer(int_kind), intent(in) :: &
1090	l, & ! level of the space-filling curve
1091	type, & ! type of SFC curve
1092	ma, & ! Major axis [0,1]
1093	md, & ! direction of major axis [-1,1]
1094	ja, & ! joiner axis [0,1]
1095	jd ! direction of joiner axis [-1,1]
1096
1097	! !OUTPUT PARAMETERS
1098
1099	integer(int_kind) :: ierr ! error return code
1100
1101	!EOP
1102	!BOC
1103	!-----------------------------------------------------------------------
1104
1105	!-------------------------------------------------
1106	! create the space-filling curve on the next level
1107	!-------------------------------------------------
1108	if(type == 2) then
1109	ierr = Hilbert(l,type,ma,md,ja,jd)
1110	elseif ( type == 3) then
1111	ierr = PeanoM(l,type,ma,md,ja,jd)
1112	elseif ( type == 5) then
1113	ierr = Cinco(l,type,ma,md,ja,jd)
1114	endif
1115
1116	!EOP
1117	!-----------------------------------------------------------------------
1118
1119	end function GenCurve
1120
1121	function FirstFactor(fac) result(res)
1122	type (factor_t) :: fac
1123	integer :: res
1124	logical :: found
1125	integer (int_kind) :: i
1126
1127	found = .false.
1128	i=1
1129	do while (i<=fac%numfact .and. (.not. found))
1130	if(fac%used(i) == 0) then
1131	res = fac%factors(i)
1132	found = .true.
1133	endif
1134	i=i+1
1135	enddo
1136
1137	end function FirstFactor
1138
1139	function FindandMark(fac,val,f2) result(found)
1140	type (factor_t) :: fac
1141	integer :: val
1142	logical :: found
1143	logical :: f2
1144	integer (int_kind) :: i
1145
1146	found = .false.
1147	i=1
1148	do while (i<=fac%numfact .and. (.not. found))
1149	if(fac%used(i) == 0) then
1150	if(fac%factors(i) .eq. val) then
1151	if(f2) then
1152	fac%used(i) = 1
1153	found = .true.
1154	else if( .not. f2) then
1155	fac%used(i) = -1
1156	found = .false.
1157	endif
1158	endif
1159	endif
1160	i=i+1
1161	enddo
1162
1163	end function FindandMark
1164
1165
1166	subroutine MatchFactor(fac1,fac2,val,found)
1167	type (factor_t) :: fac1
1168	type (factor_t) :: fac2
1169	integer :: val
1170	integer :: val1
1171	logical :: found
1172	logical :: tmp
1173
1174	found = .false.
1175
1176	val1 = FirstFactor(fac1)
1177	print *,'Matchfactor: found value: ',val
1178	found = FindandMark(fac2,val1,.true.)
1179	tmp = FindandMark(fac1,val1,found)
1180	if (found) then
1181	val = val1
1182	else
1183	val = 1
1184	endif
1185
1186	end subroutine MatchFactor
1187
1188	function ProdFactor(fac) result(res)
1189
1190	type (factor_t) :: fac
1191	integer :: res
1192	integer (int_kind) :: i
1193
1194	res = 1
1195	do i=1,fac%numfact
1196	if(fac%used(i) <= 0) then
1197	res = res * fac%factors(i)
1198	endif
1199	enddo
1200
1201	end function ProdFactor
1202
1203
1204	subroutine PrintFactor(msg,fac)
1205
1206
1207	character(len=*) :: msg
1208	type (factor_t) :: fac
1209	integer (int_kind) :: i
1210
1211	write(,) ' '
1212	write(,) 'PrintFactor: ',msg
1213	write(,) (fac%factors(i),i=1,fac%numfact)
1214	write(,) (fac%used(i),i=1,fac%numfact)
1215
1216
1217	end subroutine PrintFactor
1218
1219	!***********************************************************************
1220	!BOP
1221	! !IROUTINE: Factor
1222	! !INTERFACE:
1223
1224	function Factor(num) result(res)
1225
1226	! !DESCRIPTION:
1227	! This function factors the input value num into a
1228	! product of 2,3, and 5.
1229	!
1230	! !REVISION HISTORY:
1231	! same as module
1232	!
1233
1234	! !INPUT PARAMETERS:
1235
1236	integer(int_kind), intent(in) :: num ! number to factor
1237
1238	! !OUTPUT PARAMETERS:
1239
1240	type (factor_t) :: res
1241
1242	!EOP
1243	!BOC
1244	!-----------------------------------------------------------------------
1245	!
1246	! local variables
1247	!
1248	!-----------------------------------------------------------------------
1249
1250	integer(int_kind) :: &
1251	tmp,tmp2,tmp3,tmp5 ! tempories for the factorization algorithm
1252	integer(int_kind) :: i,n ! loop tempories
1253	logical :: found ! logical temporary
1254
1255	! --------------------------------------
1256	! Allocate allocate for max # of factors
1257	! --------------------------------------
1258	tmp = num
1259	tmp2 = log2(num)
1260	allocate(res%factors(tmp2))
1261	allocate(res%used(tmp2))
1262
1263	res%used = 0
1264	n=0
1265
1266	!-----------------------
1267	! Look for factors of 5
1268	!-----------------------
1269	found=.TRUE.
1270	do while (found)
1271	found = .FALSE.
1272	tmp5 = tmp/5
1273	if( tmp5*5 == tmp ) then
1274	n = n + 1
1275	res%factors(n) = 5
1276	found = .TRUE.
1277	tmp = tmp5
1278	endif
1279	enddo
1280
1281	!-----------------------
1282	! Look for factors of 3
1283	!-----------------------
1284	found=.TRUE.
1285	do while (found)
1286	found = .FALSE.
1287	tmp3 = tmp/3
1288	if( tmp3*3 == tmp ) then
1289	n = n + 1
1290	res%factors(n) = 3
1291	found = .TRUE.
1292	tmp = tmp3
1293	endif
1294	enddo
1295
1296	!-----------------------
1297	! Look for factors of 2
1298	!-----------------------
1299	found=.TRUE.
1300	do while (found)
1301	found = .FALSE.
1302	tmp2 = tmp/2
1303	if( tmp2*2 == tmp ) then
1304	n = n + 1
1305	res%factors(n) = 2
1306	found = .TRUE.
1307	tmp = tmp2
1308	endif
1309	enddo
1310
1311	!------------------------------------
1312	! make sure that the input value
1313	! only contains factors of 2,3,and 5
1314	!------------------------------------
1315	tmp=1
1316	do i=1,n
1317	tmp = tmp * res%factors(i)
1318	enddo
1319	if(tmp == num) then
1320	res%numfact = n
1321	else
1322	res%numfact = -1
1323	endif
1324
1325	!EOP
1326	!---------------------------------------------------------
1327	end function Factor
1328
1329	!***********************************************************************
1330	!BOP
1331	! !IROUTINE: IsFactorable
1332	! !INTERFACE:
1333
1334	function IsFactorable(n)
1335
1336	! !DESCRIPTION:
1337	! This function determines if we can factor
1338	! n into 2,3,and 5.
1339	!
1340	! !REVISION HISTORY:
1341	! same as module
1342
1343
1344	! !INTPUT PARAMETERS:
1345
1346	integer(int_kind), intent(in) :: n ! number to factor
1347
1348	! !OUTPUT PARAMETERS:
1349	logical :: IsFactorable ! .TRUE. if it is factorable
1350
1351	!EOP
1352	!BOC
1353	!-----------------------------------------------------------------------
1354	!
1355	! local variables
1356	!
1357	!-----------------------------------------------------------------------
1358
1359	type (factor_t) :: fact ! data structure to store factor information
1360
1361	fact = Factor(n)
1362	if(fact%numfact .ne. -1) then
1363	IsFactorable = .TRUE.
1364	else
1365	IsFactorable = .FALSE.
1366	endif
1367
1368	!EOP
1369	!-----------------------------------------------------------------------
1370
1371	end function IsFactorable
1372
1373	!***********************************************************************
1374	!BOP
1375	! !IROUTINE: map
1376	! !INTERFACE:
1377
1378	subroutine map(l)
1379
1380	! !DESCRIPTION:
1381	! Interface routine between internal subroutines and public
1382	! subroutines.
1383	!
1384	! !REVISION HISTORY:
1385	! same as module
1386
1387	! !INPUT PARAMETERS:
1388	integer(int_kind) :: l ! level of space-filling curve
1389
1390
1391	!EOP
1392	!BOC
1393	!-----------------------------------------------------------------------
1394	!
1395	! local variables
1396	!
1397	!-----------------------------------------------------------------------
1398
1399	integer(int_kind) :: &
1400	d, & ! dimension of curve only 2D is supported
1401	type, & ! type of space-filling curve to start off
1402	ierr ! error return code
1403
1404	d = SIZE(pos)
1405
1406	pos=0
1407	maxdim=d
1408	vcnt=0
1409
1410	type = fact%factors(l)
1411	ierr = GenCurve(l,type,0,1,0,1)
1412
1413
1414	!EOP
1415	!-----------------------------------------------------------------------
1416
1417	end subroutine map
1418
1419	!***********************************************************************
1420	!BOP
1421	! !IROUTINE: PrintCurve
1422	! !INTERFACE:
1423
1424	subroutine PrintCurve(Mesh)
1425
1426
1427	! !DESCRIPTION:
1428	! This subroutine prints the several low order
1429	! space-filling curves in an easy to read format
1430	!
1431	! !REVISION HISTORY:
1432	! same as module
1433	!
1434	! !INPUT PARAMETERS:
1435
1436	integer(int_kind), intent(in), target :: Mesh(:,:) ! SFC to be printed
1437
1438	!EOP
1439	!BOC
1440	!-----------------------------------------------------------------------
1441	!
1442	! local variables
1443	!
1444	!-----------------------------------------------------------------------
1445	integer(int_kind) :: &
1446	gridsize, &! order of space-filling curve
1447	i ! loop temporary
1448
1449	!-----------------------------------------------------------------------
1450
1451	gridsize = SIZE(Mesh,dim=1)
1452
1453	if(gridsize == 2) then
1454	write (,) "A Level 1 Hilbert Curve:"
1455	write (,) "------------------------"
1456	do i=1,gridsize
1457	write(*,2) Mesh(1,i),Mesh(2,i)
1458	enddo
1459	else if(gridsize == 3) then
1460	write (,) "A Level 1 Peano Meandering Curve:"
1461	write (,) "---------------------------------"
1462	do i=1,gridsize
1463	write(*,3) Mesh(1,i),Mesh(2,i),Mesh(3,i)
1464	enddo
1465	else if(gridsize == 4) then
1466	write (,) "A Level 2 Hilbert Curve:"
1467	write (,) "------------------------"
1468	do i=1,gridsize
1469	write(*,4) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i)
1470	enddo
1471	else if(gridsize == 5) then
1472	write (,) "A Level 1 Cinco Curve:"
1473	write (,) "------------------------"
1474	do i=1,gridsize
1475	write(*,5) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i),Mesh(5,i)
1476	enddo
1477	else if(gridsize == 6) then
1478	write (,) "A Level 1 Hilbert and Level 1 Peano Curve:"
1479	write (,) "------------------------------------------"
1480	do i=1,gridsize
1481	write(*,6) Mesh(1,i),Mesh(2,i),Mesh(3,i), &
1482	Mesh(4,i),Mesh(5,i),Mesh(6,i)
1483	enddo
1484	else if(gridsize == 8) then
1485	write (,) "A Level 3 Hilbert Curve:"
1486	write (,) "------------------------"
1487	do i=1,gridsize
1488	write(*,8) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1489	Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i)
1490	enddo
1491	else if(gridsize == 9) then
1492	write (,) "A Level 2 Peano Meandering Curve:"
1493	write (,) "---------------------------------"
1494	do i=1,gridsize
1495	write(*,9) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1496	Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1497	Mesh(9,i)
1498	enddo
1499	else if(gridsize == 10) then
1500	write (,) "A Level 1 Hilbert and Level 1 Cinco Curve:"
1501	write (,) "---------------------------------"
1502	do i=1,gridsize
1503	write(*,10) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1504	Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1505	Mesh(9,i),Mesh(10,i)
1506	enddo
1507	else if(gridsize == 12) then
1508	write (,) "A Level 2 Hilbert and Level 1 Peano Curve:"
1509	write (,) "------------------------------------------"
1510	do i=1,gridsize
1511	write(*,12) Mesh(1,i),Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1512	Mesh(5,i),Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1513	Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i)
1514	enddo
1515	else if(gridsize == 15) then
1516	write (,) "A Level 1 Peano and Level 1 Cinco Curve:"
1517	write (,) "------------------------"
1518	do i=1,gridsize
1519	write(*,15) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1520	Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1521	Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1522	Mesh(13,i),Mesh(14,i),Mesh(15,i)
1523	enddo
1524	else if(gridsize == 16) then
1525	write (,) "A Level 4 Hilbert Curve:"
1526	write (,) "------------------------"
1527	do i=1,gridsize
1528	write(*,16) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1529	Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1530	Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1531	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i)
1532	enddo
1533	else if(gridsize == 18) then
1534	write (,) "A Level 1 Hilbert and Level 2 Peano Curve:"
1535	write (,) "------------------------------------------"
1536	do i=1,gridsize
1537	write(*,18) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1538	Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1539	Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1540	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1541	Mesh(17,i),Mesh(18,i)
1542	enddo
1543	else if(gridsize == 20) then
1544	write (,) "A Level 2 Hilbert and Level 1 Cinco Curve:"
1545	write (,) "------------------------------------------"
1546	do i=1,gridsize
1547	write(*,20) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1548	Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1549	Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1550	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1551	Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i)
1552	enddo
1553	else if(gridsize == 24) then
1554	write (,) "A Level 3 Hilbert and Level 1 Peano Curve:"
1555	write (,) "------------------------------------------"
1556	do i=1,gridsize
1557	write(*,24) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1558	Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1559	Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1560	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1561	Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1562	Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i)
1563	enddo
1564	else if(gridsize == 25) then
1565	write (,) "A Level 2 Cinco Curve:"
1566	write (,) "------------------------------------------"
1567	do i=1,gridsize
1568	write(*,25) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1569	Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1570	Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1571	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1572	Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1573	Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
1574	Mesh(25,i)
1575	enddo
1576	else if(gridsize == 27) then
1577	write (,) "A Level 3 Peano Meandering Curve:"
1578	write (,) "---------------------------------"
1579	do i=1,gridsize
1580	write(*,27) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1581	Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1582	Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1583	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1584	Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1585	Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
1586	Mesh(25,i),Mesh(26,i),Mesh(27,i)
1587	enddo
1588	else if(gridsize == 32) then
1589	write (,) "A Level 5 Hilbert Curve:"
1590	write (,) "------------------------"
1591	do i=1,gridsize
1592	write(*,32) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1593	Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1594	Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1595	Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1596	Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1597	Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
1598	Mesh(25,i),Mesh(26,i),Mesh(27,i),Mesh(28,i), &
1599	Mesh(29,i),Mesh(30,i),Mesh(31,i),Mesh(32,i)
1600	enddo
1601	endif
1602	2 format('\|',2(i2,'\|'))
1603	3 format('\|',3(i2,'\|'))
1604	4 format('\|',4(i2,'\|'))
1605	5 format('\|',5(i2,'\|'))
1606	6 format('\|',6(i2,'\|'))
1607	8 format('\|',8(i2,'\|'))
1608	9 format('\|',9(i2,'\|'))
1609	10 format('\|',10(i2,'\|'))
1610	12 format('\|',12(i3,'\|'))
1611	15 format('\|',15(i3,'\|'))
1612	16 format('\|',16(i3,'\|'))
1613	18 format('\|',18(i3,'\|'))
1614	20 format('\|',20(i3,'\|'))
1615	24 format('\|',24(i3,'\|'))
1616	25 format('\|',25(i3,'\|'))
1617	27 format('\|',27(i3,'\|'))
1618	32 format('\|',32(i4,'\|'))
1619
1620	!EOC
1621	!-----------------------------------------------------------------------
1622
1623	end subroutine PrintCurve
1624
1625	!***********************************************************************
1626	!BOP
1627	! !IROUTINE: GenSpaceCurve
1628	! !INTERFACE:
1629
1630	subroutine GenSpaceCurve(Mesh)
1631
1632	! !DESCRIPTION:
1633	! This subroutine is the public interface into the
1634	! space-filling curve functionality
1635	!
1636	! !REVISION HISTORY:
1637	! same as module
1638	!
1639
1640	! !INPUT/OUTPUT PARAMETERS:
1641	integer(int_kind), target,intent(inout) :: &
1642	Mesh(:,:) ! The SFC ordering in 2D array
1643
1644	#if 0
1645	integer(int_kind), intent(inout) :: &
1646	Path(:,:) ! The path for the domain traversal
1647	#endif
1648
1649	!EOP
1650	!BOC
1651	!-----------------------------------------------------------------------
1652	!
1653	! local variables
1654	!
1655	!-----------------------------------------------------------------------
1656
1657	integer(int_kind) :: &
1658	level, &! Level of space-filling curve
1659	dim ! dimension of SFC... currently limited to 2D
1660
1661	integer(int_kind) :: gridsize ! number of points on a side
1662
1663	!-----------------------------------------------------------------------
1664
1665	!-----------------------------------------
1666	! Setup the size of the grid to traverse
1667	!-----------------------------------------
1668	dim = 2
1669	gridsize = SIZE(Mesh,dim=1)
1670	fact = factor(gridsize)
1671	level = fact%numfact
1672
1673	if(verbose) write(,) 'GenSpacecurve: level is ',level
1674	allocate(ordered(gridsize,gridsize))
1675
1676	!--------------------------------------------
1677	! Setup the working arrays for the traversal
1678	!--------------------------------------------
1679	allocate(pos(0:dim-1))
1680
1681	!-----------------------------------------------------
1682	! The array ordered will contain the visitation order
1683	!-----------------------------------------------------
1684	ordered(:,:) = 0
1685
1686	call map(level)
1687
1688	Mesh(:,:) = ordered(:,:)
1689
1690	deallocate(pos,ordered)
1691
1692	!EOP
1693	!-----------------------------------------------------------------------
1694
1695	end subroutine GenSpaceCurve
1696
1697	end module spacecurve_mod
1698
1699	!\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|\|

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format