1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
/***************************************************************************
* blitz/traversal.cc Space-filling curve based traversal orders
*
* $Id: traversal.cc 1414 2005-11-01 22:04:59Z cookedm $
*
* Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org>
*
* This code was relicensed under the modified BSD license for use in SciPy
* by Todd Veldhuizen (see LICENSE.txt in the weave directory).
*
*
* Suggestions: blitz-dev@oonumerics.org
* Bugs: blitz-bugs@oonumerics.org
*
* For more information, please see the Blitz++ Home Page:
* http://oonumerics.org/blitz/
*
***************************************************************************/
#ifndef BZ_TRAVERSAL_CC
#define BZ_TRAVERSAL_CC
#ifndef BZ_TRAVERSAL_H
#error <blitz/traversal.cc> must be included via <blitz/traversal.h>
#endif
BZ_NAMESPACE(blitz)
// Next line is a workaround for Intel C++ V4.0 oddity, due
// to Allan Stokes.
static set<TraversalOrder<2> > *_bz_intel_kludge;
//template<int N_dimensions>
//_bz_typename TraversalOrderCollection<N_dimensions>::T_set
// TraversalOrderCollection<N_dimensions>::traversals_;
template<int N>
void makeHilbert(Vector<TinyVector<int,N> >& coord,
int x0, int y0, int xis, int xjs,
int yis, int yjs, int n, int& i)
{
// N != 2 is not yet implemented.
BZPRECONDITION(N == 2);
if (!n)
{
if (i > coord.length())
{
cerr << "makeHilbert: vector not long enough" << endl;
exit(1);
}
coord[i][0] = x0 + (xis+yis)/2;
coord[i][1] = y0 + (xjs+yjs)/2;
++i;
}
else {
makeHilbert(coord,x0,y0,yis/2, yjs/2, xis/2, xjs/2, n-1, i);
makeHilbert(coord,x0+xis/2,y0+xjs/2,xis/2,xjs/2,yis/2,yjs/2,n-1,i);
makeHilbert(coord,x0+xis/2+yis/2,y0+xjs/2+yjs/2,xis/2,xjs/2,yis/2,
yjs/2,n-1,i);
makeHilbert(coord,x0+xis/2+yis, y0+xjs/2+yjs, -yis/2,-yjs/2,-xis/2,
-xjs/2,n-1,i);
}
}
template<int N_dimensions>
void MakeHilbertTraversal(Vector<TinyVector<int,N_dimensions> >& coord,
int length)
{
// N_dimensions != 2 not yet implemented
BZPRECONDITION(N_dimensions == 2);
// The constant on the next line is ln(2)
int d = (int)(::ceil(::log((double)length) / 0.693147180559945309417));
int N = 1 << d;
const int Npoints = N*N;
Vector<TinyVector<int,2> > coord2(Npoints);
int i=0;
makeHilbert(coord2,0,0,32768,0,0,32768,d,i);
int xp0 = coord2[0][0];
int yp0 = coord2[0][1];
int j=0;
coord.resize(length * length);
for (int i=0; i < Npoints; ++i)
{
coord2[i][0] = (coord2[i][0]-xp0)/(2*xp0);
coord2[i][1] = (coord2[i][1]-yp0)/(2*yp0);
if ((coord2[i][0] < length) && (coord2[i][1] < length)
&& (coord2[i][0] >= 0) && (coord2[i][1] >= 0))
{
coord[j][0] = coord2[i][0];
coord[j][1] = coord2[i][1];
++j;
}
}
}
template<int N_dimensions>
void generateFastTraversalOrder(const TinyVector<int,N_dimensions>& size)
{
BZPRECONDITION(N_dimensions == 2);
BZPRECONDITION(size[0] == size[1]);
TraversalOrderCollection<2> travCol;
if (travCol.find(size))
return;
Vector<TinyVector<int,2> > ordering(size[0]);
MakeHilbertTraversal(ordering, size[0]);
travCol.insert(TraversalOrder<2>(size, ordering));
}
BZ_NAMESPACE_END
#endif // BZ_TRAVERSAL_CC
|
