Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88394994
peano.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 18, 14:14
Size
12 KB
Mime Type
text/x-c
Expires
Sun, Oct 20, 14:14 (2 d)
Engine
blob
Format
Raw Data
Handle
21764794
Attached To
R195 SCITAS application benchmarks
peano.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file peano.c
* \brief Routines to compute a Peano-Hilbert order
*
* This file contains routines to compute Peano-Hilbert keys, and to put the
* particle data into the order of these keys, i.e. into the order of a
* space-filling fractal curve.
*/
static
struct
peano_hilbert_data
{
peanokey
key
;
int
index
;
}
*
mp
;
static
int
*
Id
;
/*! This function puts the particles into Peano-Hilbert order by sorting them
* according to their keys. The latter half already been computed in the
* domain decomposition. Since gas particles need to stay at the beginning of
* the particle list, they are sorted as a separate block.
*/
void
peano_hilbert_order
(
void
)
{
int
i
;
if
(
ThisTask
==
0
)
printf
(
"begin Peano-Hilbert order...
\n
"
);
if
(
N_gas
)
{
mp
=
malloc
(
sizeof
(
struct
peano_hilbert_data
)
*
N_gas
);
Id
=
malloc
(
sizeof
(
int
)
*
N_gas
);
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
mp
[
i
].
index
=
i
;
mp
[
i
].
key
=
Key
[
i
];
}
qsort
(
mp
,
N_gas
,
sizeof
(
struct
peano_hilbert_data
),
compare_key
);
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
Id
[
mp
[
i
].
index
]
=
i
;
reorder_gas
();
free
(
Id
);
free
(
mp
);
}
#ifdef STELLAR_PROP
if
(
N_stars
>
0
)
{
mp
=
malloc
(
sizeof
(
struct
peano_hilbert_data
)
*
(
N_stars
));
mp
-=
(
N_gas
);
Id
=
malloc
(
sizeof
(
int
)
*
(
N_stars
));
Id
-=
(
N_gas
);
for
(
i
=
N_gas
;
i
<
N_gas
+
N_stars
;
i
++
)
{
mp
[
i
].
index
=
i
;
mp
[
i
].
key
=
Key
[
i
];
}
qsort
(
mp
+
N_gas
,
N_stars
,
sizeof
(
struct
peano_hilbert_data
),
compare_key
);
for
(
i
=
N_gas
;
i
<
N_gas
+
N_stars
;
i
++
)
Id
[
mp
[
i
].
index
]
=
i
;
reorder_stars
();
Id
+=
N_gas
;
free
(
Id
);
mp
+=
N_gas
;
free
(
mp
);
}
if
(
NumPart
-
N_gas
-
N_stars
>
0
)
{
mp
=
malloc
(
sizeof
(
struct
peano_hilbert_data
)
*
(
NumPart
-
N_gas
-
N_stars
));
mp
-=
(
N_gas
+
N_stars
);
Id
=
malloc
(
sizeof
(
int
)
*
(
NumPart
-
N_gas
-
N_stars
));
Id
-=
(
N_gas
+
N_stars
);
for
(
i
=
N_gas
+
N_stars
;
i
<
NumPart
;
i
++
)
{
mp
[
i
].
index
=
i
;
mp
[
i
].
key
=
Key
[
i
];
}
qsort
(
mp
+
N_gas
+
N_stars
,
NumPart
-
N_gas
-
N_stars
,
sizeof
(
struct
peano_hilbert_data
),
compare_key
);
for
(
i
=
N_gas
+
N_stars
;
i
<
NumPart
;
i
++
)
Id
[
mp
[
i
].
index
]
=
i
;
reorder_particles
();
Id
+=
N_gas
+
N_stars
;
free
(
Id
);
mp
+=
N_gas
+
N_stars
;
free
(
mp
);
}
/* now, do StP */
if
(
N_stars
>
0
)
{
Id
=
malloc
(
sizeof
(
int
)
*
(
N_stars
));
for
(
i
=
N_gas
;
i
<
N_gas
+
N_stars
;
i
++
)
{
Id
[
P
[
i
].
StPIdx
]
=
i
-
N_gas
;
}
reorder_st
();
free
(
Id
);
}
#else
if
(
NumPart
-
N_gas
>
0
)
{
mp
=
malloc
(
sizeof
(
struct
peano_hilbert_data
)
*
(
NumPart
-
N_gas
));
mp
-=
(
N_gas
);
Id
=
malloc
(
sizeof
(
int
)
*
(
NumPart
-
N_gas
));
Id
-=
(
N_gas
);
for
(
i
=
N_gas
;
i
<
NumPart
;
i
++
)
{
mp
[
i
].
index
=
i
;
mp
[
i
].
key
=
Key
[
i
];
}
qsort
(
mp
+
N_gas
,
NumPart
-
N_gas
,
sizeof
(
struct
peano_hilbert_data
),
compare_key
);
for
(
i
=
N_gas
;
i
<
NumPart
;
i
++
)
Id
[
mp
[
i
].
index
]
=
i
;
reorder_particles
();
Id
+=
N_gas
;
free
(
Id
);
mp
+=
N_gas
;
free
(
mp
);
}
#endif
#ifdef CHECK_ID_CORRESPONDENCE
if
(
ThisTask
==
0
)
printf
(
"Check id correspondence after peano-hilbert order...
\n
"
);
for
(
i
=
N_gas
;
i
<
N_gas
+
N_stars
;
i
++
)
{
//printf("-> i=%d ID=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n",i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
if
(
P
[
i
].
ID
!=
StP
[
i
-
N_gas
].
ID
)
{
printf
(
"
\n
P/StP ID correspondance error
\n
"
);
printf
(
"(%d) (in peano-hilbert) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d StP[P[i].StPIdx].ID=%d
\n\n
"
,
ThisTask
,
N_stars
,
N_gas
,
i
,
P
[
i
].
ID
,
P
[
i
].
StPIdx
,
StP
[
P
[
i
].
StPIdx
].
PIdx
,
StP
[
P
[
i
].
StPIdx
].
ID
);
endrun
(
1212001
);
}
if
(
StP
[
P
[
i
].
StPIdx
].
PIdx
!=
i
)
{
printf
(
"
\n
P/StP correspondance error
\n
"
);
printf
(
"(%d) (in peano-hilbert) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d
\n\n
"
,
ThisTask
,
N_stars
,
N_gas
,
i
,
P
[
i
].
ID
,
P
[
i
].
StPIdx
,
StP
[
P
[
i
].
StPIdx
].
PIdx
);
endrun
(
1212001
);
}
if
(
StP
[
P
[
i
].
StPIdx
].
ID
!=
P
[
i
].
ID
)
{
printf
(
"
\n
P/StP correspondance error
\n
"
);
printf
(
"(%d) (in peano-hilbert) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d
\n\n
"
,
ThisTask
,
N_gas
,
N_stars
,
i
,
P
[
i
].
Type
,
P
[
i
].
ID
,
P
[
i
].
StPIdx
,
StP
[
P
[
i
].
StPIdx
].
ID
);
endrun
(
1212002
);
}
}
if
(
ThisTask
==
0
)
printf
(
"Check id correspondence after peano-hilbert order...
\n
"
);
#endif
if
(
ThisTask
==
0
)
printf
(
"Peano-Hilbert done.
\n
"
);
}
/*! This function is a comparison kernel for sorting the Peano-Hilbert keys.
*/
int
compare_key
(
const
void
*
a
,
const
void
*
b
)
{
if
(((
struct
peano_hilbert_data
*
)
a
)
->
key
<
(((
struct
peano_hilbert_data
*
)
b
)
->
key
))
return
-
1
;
if
(((
struct
peano_hilbert_data
*
)
a
)
->
key
>
(((
struct
peano_hilbert_data
*
)
b
)
->
key
))
return
+
1
;
return
0
;
}
/*! This function brings the gas particles into the same order as the sorted
* keys. (The sort is first done only on the keys themselves and done
* directly on the gas particles in order to reduce the amount of data that
* needs to be moved in memory. Only once the order is established, the gas
* particles are rearranged, such that each particle has to be moved at most
* once.)
*/
void
reorder_gas
(
void
)
{
int
i
;
struct
particle_data
Psave
,
Psource
;
struct
sph_particle_data
SphPsave
,
SphPsource
;
int
idsource
,
idsave
,
dest
;
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
if
(
Id
[
i
]
!=
i
)
{
Psource
=
P
[
i
];
SphPsource
=
SphP
[
i
];
idsource
=
Id
[
i
];
dest
=
Id
[
i
];
do
{
Psave
=
P
[
dest
];
SphPsave
=
SphP
[
dest
];
idsave
=
Id
[
dest
];
P
[
dest
]
=
Psource
;
SphP
[
dest
]
=
SphPsource
;
Id
[
dest
]
=
idsource
;
if
(
dest
==
i
)
break
;
Psource
=
Psave
;
SphPsource
=
SphPsave
;
idsource
=
idsave
;
dest
=
idsource
;
}
while
(
1
);
}
}
}
#ifdef STELLAR_PROP
/*! This function brings the stars into the same order as
* the sorted keys. (The sort is first done only on the keys themselves and
* done directly on the particles in order to reduce the amount of data that
* needs to be moved in memory. Only once the order is established, the
* particles are rearranged, such that each particle has to be moved at most
* once.)
*/
void
reorder_stars
(
void
)
{
int
i
;
struct
particle_data
Psave
,
Psource
;
int
idsource
,
idsave
,
dest
;
for
(
i
=
N_gas
;
i
<
N_gas
+
N_stars
;
i
++
)
{
if
(
Id
[
i
]
!=
i
)
{
Psource
=
P
[
i
];
idsource
=
Id
[
i
];
dest
=
Id
[
i
];
do
{
Psave
=
P
[
dest
];
idsave
=
Id
[
dest
];
P
[
dest
]
=
Psource
;
Id
[
dest
]
=
idsource
;
/* restore the link with Stp */
StP
[
P
[
dest
].
StPIdx
].
PIdx
=
dest
;
if
(
dest
==
i
)
break
;
Psource
=
Psave
;
idsource
=
idsave
;
dest
=
idsource
;
}
while
(
1
);
}
}
}
void
reorder_st
(
void
)
{
int
i
;
struct
st_particle_data
StPsave
,
StPsource
;
int
idsource
,
idsave
,
dest
;
for
(
i
=
0
;
i
<
N_stars
;
i
++
)
{
if
(
Id
[
i
]
!=
i
)
{
StPsource
=
StP
[
i
];
idsource
=
Id
[
i
];
dest
=
Id
[
i
];
do
{
StPsave
=
StP
[
dest
];
idsave
=
Id
[
dest
];
StP
[
dest
]
=
StPsource
;
Id
[
dest
]
=
idsource
;
/* restore the link with P */
P
[
StP
[
dest
].
PIdx
].
StPIdx
=
dest
;
if
(
dest
==
i
)
break
;
StPsource
=
StPsave
;
idsource
=
idsave
;
dest
=
idsource
;
}
while
(
1
);
}
}
}
#endif
/*! This function brings the collisionless particles into the same order as
* the sorted keys. (The sort is first done only on the keys themselves and
* done directly on the particles in order to reduce the amount of data that
* needs to be moved in memory. Only once the order is established, the
* particles are rearranged, such that each particle has to be moved at most
* once.)
*/
void
reorder_particles
(
void
)
{
int
i
;
struct
particle_data
Psave
,
Psource
;
int
idsource
,
idsave
,
dest
;
#ifdef STELLAR_PROP
for
(
i
=
N_gas
+
N_stars
;
i
<
NumPart
;
i
++
)
#else
for
(
i
=
N_gas
;
i
<
NumPart
;
i
++
)
#endif
{
if
(
Id
[
i
]
!=
i
)
{
Psource
=
P
[
i
];
idsource
=
Id
[
i
];
dest
=
Id
[
i
];
do
{
Psave
=
P
[
dest
];
idsave
=
Id
[
dest
];
P
[
dest
]
=
Psource
;
Id
[
dest
]
=
idsource
;
if
(
dest
==
i
)
break
;
Psource
=
Psave
;
idsource
=
idsave
;
dest
=
idsource
;
}
while
(
1
);
}
}
}
static
int
quadrants
[
24
][
2
][
2
][
2
]
=
{
/* rotx=0, roty=0-3 */
{{{
0
,
7
},
{
1
,
6
}},
{{
3
,
4
},
{
2
,
5
}}},
{{{
7
,
4
},
{
6
,
5
}},
{{
0
,
3
},
{
1
,
2
}}},
{{{
4
,
3
},
{
5
,
2
}},
{{
7
,
0
},
{
6
,
1
}}},
{{{
3
,
0
},
{
2
,
1
}},
{{
4
,
7
},
{
5
,
6
}}},
/* rotx=1, roty=0-3 */
{{{
1
,
0
},
{
6
,
7
}},
{{
2
,
3
},
{
5
,
4
}}},
{{{
0
,
3
},
{
7
,
4
}},
{{
1
,
2
},
{
6
,
5
}}},
{{{
3
,
2
},
{
4
,
5
}},
{{
0
,
1
},
{
7
,
6
}}},
{{{
2
,
1
},
{
5
,
6
}},
{{
3
,
0
},
{
4
,
7
}}},
/* rotx=2, roty=0-3 */
{{{
6
,
1
},
{
7
,
0
}},
{{
5
,
2
},
{
4
,
3
}}},
{{{
1
,
2
},
{
0
,
3
}},
{{
6
,
5
},
{
7
,
4
}}},
{{{
2
,
5
},
{
3
,
4
}},
{{
1
,
6
},
{
0
,
7
}}},
{{{
5
,
6
},
{
4
,
7
}},
{{
2
,
1
},
{
3
,
0
}}},
/* rotx=3, roty=0-3 */
{{{
7
,
6
},
{
0
,
1
}},
{{
4
,
5
},
{
3
,
2
}}},
{{{
6
,
5
},
{
1
,
2
}},
{{
7
,
4
},
{
0
,
3
}}},
{{{
5
,
4
},
{
2
,
3
}},
{{
6
,
7
},
{
1
,
0
}}},
{{{
4
,
7
},
{
3
,
0
}},
{{
5
,
6
},
{
2
,
1
}}},
/* rotx=4, roty=0-3 */
{{{
6
,
7
},
{
5
,
4
}},
{{
1
,
0
},
{
2
,
3
}}},
{{{
7
,
0
},
{
4
,
3
}},
{{
6
,
1
},
{
5
,
2
}}},
{{{
0
,
1
},
{
3
,
2
}},
{{
7
,
6
},
{
4
,
5
}}},
{{{
1
,
6
},
{
2
,
5
}},
{{
0
,
7
},
{
3
,
4
}}},
/* rotx=5, roty=0-3 */
{{{
2
,
3
},
{
1
,
0
}},
{{
5
,
4
},
{
6
,
7
}}},
{{{
3
,
4
},
{
0
,
7
}},
{{
2
,
5
},
{
1
,
6
}}},
{{{
4
,
5
},
{
7
,
6
}},
{{
3
,
2
},
{
0
,
1
}}},
{{{
5
,
2
},
{
6
,
1
}},
{{
4
,
3
},
{
7
,
0
}}}
};
static
int
rotxmap_table
[
24
]
=
{
4
,
5
,
6
,
7
,
8
,
9
,
10
,
11
,
12
,
13
,
14
,
15
,
0
,
1
,
2
,
3
,
17
,
18
,
19
,
16
,
23
,
20
,
21
,
22
};
static
int
rotymap_table
[
24
]
=
{
1
,
2
,
3
,
0
,
16
,
17
,
18
,
19
,
11
,
8
,
9
,
10
,
22
,
23
,
20
,
21
,
14
,
15
,
12
,
13
,
4
,
5
,
6
,
7
};
static
int
rotx_table
[
8
]
=
{
3
,
0
,
0
,
2
,
2
,
0
,
0
,
1
};
static
int
roty_table
[
8
]
=
{
0
,
1
,
1
,
2
,
2
,
3
,
3
,
0
};
static
int
sense_table
[
8
]
=
{
-
1
,
-
1
,
-
1
,
+
1
,
+
1
,
-
1
,
-
1
,
-
1
};
static
int
flag_quadrants_inverse
=
1
;
static
char
quadrants_inverse_x
[
24
][
8
];
static
char
quadrants_inverse_y
[
24
][
8
];
static
char
quadrants_inverse_z
[
24
][
8
];
/*! This function computes a Peano-Hilbert key for an integer triplet (x,y,z),
* with x,y,z in the range between 0 and 2^bits-1.
*/
peanokey
peano_hilbert_key
(
int
x
,
int
y
,
int
z
,
int
bits
)
{
int
i
,
quad
,
bitx
,
bity
,
bitz
;
int
mask
,
rotation
,
rotx
,
roty
,
sense
;
peanokey
key
;
mask
=
1
<<
(
bits
-
1
);
key
=
0
;
rotation
=
0
;
sense
=
1
;
for
(
i
=
0
;
i
<
bits
;
i
++
,
mask
>>=
1
)
{
bitx
=
(
x
&
mask
)
?
1
:
0
;
bity
=
(
y
&
mask
)
?
1
:
0
;
bitz
=
(
z
&
mask
)
?
1
:
0
;
quad
=
quadrants
[
rotation
][
bitx
][
bity
][
bitz
];
key
<<=
3
;
key
+=
(
sense
==
1
)
?
(
quad
)
:
(
7
-
quad
);
rotx
=
rotx_table
[
quad
];
roty
=
roty_table
[
quad
];
sense
*=
sense_table
[
quad
];
while
(
rotx
>
0
)
{
rotation
=
rotxmap_table
[
rotation
];
rotx
--
;
}
while
(
roty
>
0
)
{
rotation
=
rotymap_table
[
rotation
];
roty
--
;
}
}
return
key
;
}
/*! This function computes for a given Peano-Hilbert key, the inverse,
* i.e. the integer triplet (x,y,z) with a Peano-Hilbert key equal to the
* input key. (This functionality is actually not needed in the present
* code.)
*/
void
peano_hilbert_key_inverse
(
peanokey
key
,
int
bits
,
int
*
x
,
int
*
y
,
int
*
z
)
{
int
i
,
keypart
,
bitx
,
bity
,
bitz
,
mask
,
quad
,
rotation
,
shift
;
char
sense
,
rotx
,
roty
;
if
(
flag_quadrants_inverse
)
{
flag_quadrants_inverse
=
0
;
for
(
rotation
=
0
;
rotation
<
24
;
rotation
++
)
for
(
bitx
=
0
;
bitx
<
2
;
bitx
++
)
for
(
bity
=
0
;
bity
<
2
;
bity
++
)
for
(
bitz
=
0
;
bitz
<
2
;
bitz
++
)
{
quad
=
quadrants
[
rotation
][
bitx
][
bity
][
bitz
];
quadrants_inverse_x
[
rotation
][
quad
]
=
bitx
;
quadrants_inverse_y
[
rotation
][
quad
]
=
bity
;
quadrants_inverse_z
[
rotation
][
quad
]
=
bitz
;
}
}
shift
=
3
*
(
bits
-
1
);
mask
=
7
<<
shift
;
rotation
=
0
;
sense
=
1
;
*
x
=
*
y
=
*
z
=
0
;
for
(
i
=
0
;
i
<
bits
;
i
++
,
mask
>>=
3
,
shift
-=
3
)
{
keypart
=
(
key
&
mask
)
>>
shift
;
quad
=
(
sense
==
1
)
?
(
keypart
)
:
(
7
-
keypart
);
*
x
=
(
*
x
<<
1
)
+
quadrants_inverse_x
[
rotation
][
quad
];
*
y
=
(
*
y
<<
1
)
+
quadrants_inverse_y
[
rotation
][
quad
];
*
z
=
(
*
z
<<
1
)
+
quadrants_inverse_z
[
rotation
][
quad
];
rotx
=
rotx_table
[
quad
];
roty
=
roty_table
[
quad
];
sense
*=
sense_table
[
quad
];
while
(
rotx
>
0
)
{
rotation
=
rotxmap_table
[
rotation
];
rotx
--
;
}
while
(
roty
>
0
)
{
rotation
=
rotymap_table
[
rotation
];
roty
--
;
}
}
}
Event Timeline
Log In to Comment