1 |
dl |
1.1 |
/* |
2 |
|
|
* Written by Doug Lea with assistance from members of JCP JSR-166 |
3 |
|
|
* Expert Group and released to the public domain, as explained at |
4 |
|
|
* http://creativecommons.org/licenses/publicdomain |
5 |
|
|
*/ |
6 |
|
|
|
7 |
|
|
// Adapted from a cilk benchmark |
8 |
|
|
|
9 |
|
|
import java.util.concurrent.*; |
10 |
|
|
|
11 |
|
|
public class Heat { |
12 |
|
|
static final long NPS = (1000L * 1000 * 1000); |
13 |
|
|
|
14 |
|
|
// Parameters |
15 |
|
|
static int nx; |
16 |
|
|
static int ny; |
17 |
|
|
static int nt; |
18 |
|
|
static int leafmaxcol; |
19 |
|
|
|
20 |
|
|
// the matrix representing the cells |
21 |
|
|
static double[][] newm; |
22 |
|
|
|
23 |
|
|
// alternating workspace matrix |
24 |
|
|
static double[][] oldm; |
25 |
|
|
|
26 |
|
|
public static void main(String[] args) throws Exception { |
27 |
|
|
int procs = 0; |
28 |
|
|
nx = 4096; |
29 |
|
|
ny = 1024; |
30 |
|
|
nt = 1000; |
31 |
|
|
leafmaxcol = 16; |
32 |
|
|
|
33 |
|
|
try { |
34 |
|
|
if (args.length > 0) |
35 |
|
|
procs = Integer.parseInt(args[0]); |
36 |
|
|
if (args.length > 1) |
37 |
|
|
nx = Integer.parseInt(args[1]); |
38 |
|
|
if (args.length > 2) |
39 |
|
|
ny = Integer.parseInt(args[2]); |
40 |
|
|
if (args.length > 3) |
41 |
|
|
nt = Integer.parseInt(args[3]); |
42 |
|
|
if (args.length > 4) |
43 |
|
|
leafmaxcol = Integer.parseInt(args[4]); |
44 |
|
|
} |
45 |
|
|
catch (Exception e) { |
46 |
|
|
System.out.println("Usage: java Heat threads rows cols steps granularity"); |
47 |
|
|
return; |
48 |
|
|
} |
49 |
|
|
|
50 |
jsr166 |
1.2 |
ForkJoinPool g = procs == 0? new ForkJoinPool() : |
51 |
dl |
1.1 |
new ForkJoinPool(procs); |
52 |
jsr166 |
1.2 |
|
53 |
dl |
1.1 |
System.out.print("parallelism = " + g.getParallelism()); |
54 |
|
|
System.out.print(" granularity = " + leafmaxcol); |
55 |
|
|
System.out.print(" rows = " + nx); |
56 |
|
|
System.out.print(" columns = " + ny); |
57 |
|
|
System.out.println(" steps = " + nt); |
58 |
jsr166 |
1.2 |
|
59 |
dl |
1.1 |
oldm = new double[nx][ny]; |
60 |
|
|
newm = new double[nx][ny]; |
61 |
|
|
|
62 |
|
|
for (int i = 0; i < 5; ++i) { |
63 |
|
|
long last = System.nanoTime(); |
64 |
|
|
RecursiveAction main = new RecursiveAction() { |
65 |
|
|
public void compute() { |
66 |
|
|
for (int timestep = 0; timestep <= nt; timestep++) { |
67 |
|
|
(new Compute(0, nx, timestep)).invoke(); |
68 |
|
|
} |
69 |
|
|
} |
70 |
|
|
}; |
71 |
|
|
g.invoke(main); |
72 |
|
|
double elapsed = elapsedTime(last); |
73 |
|
|
System.out.printf("time: %7.3f", elapsed); |
74 |
|
|
System.out.println(); |
75 |
|
|
} |
76 |
|
|
System.out.println(g); |
77 |
|
|
g.shutdown(); |
78 |
|
|
} |
79 |
|
|
|
80 |
|
|
static double elapsedTime(long startTime) { |
81 |
|
|
return (double)(System.nanoTime() - startTime) / NPS; |
82 |
|
|
} |
83 |
|
|
|
84 |
|
|
// constants (at least for this demo) |
85 |
|
|
static final double xu = 0.0; |
86 |
|
|
static final double xo = 1.570796326794896558; |
87 |
|
|
static final double yu = 0.0; |
88 |
|
|
static final double yo = 1.570796326794896558; |
89 |
|
|
static final double tu = 0.0; |
90 |
|
|
static final double to = 0.0000001; |
91 |
|
|
|
92 |
|
|
static final double dx = (xo - xu) / (nx - 1); |
93 |
|
|
static final double dy = (yo - yu) / (ny - 1); |
94 |
jsr166 |
1.2 |
static final double dt = (to - tu) / nt; |
95 |
dl |
1.1 |
static final double dtdxsq = dt / (dx * dx); |
96 |
|
|
static final double dtdysq = dt / (dy * dy); |
97 |
|
|
|
98 |
|
|
|
99 |
|
|
// the function being applied across the cells |
100 |
jsr166 |
1.2 |
static final double f(double x, double y) { |
101 |
|
|
return Math.sin(x) * Math.sin(y); |
102 |
dl |
1.1 |
} |
103 |
|
|
|
104 |
|
|
// random starting values |
105 |
|
|
|
106 |
jsr166 |
1.2 |
static final double randa(double x, double t) { |
107 |
|
|
return 0.0; |
108 |
dl |
1.1 |
} |
109 |
jsr166 |
1.2 |
static final double randb(double x, double t) { |
110 |
|
|
return Math.exp(-2*t) * Math.sin(x); |
111 |
dl |
1.1 |
} |
112 |
jsr166 |
1.2 |
static final double randc(double y, double t) { |
113 |
|
|
return 0.0; |
114 |
dl |
1.1 |
} |
115 |
jsr166 |
1.2 |
static final double randd(double y, double t) { |
116 |
|
|
return Math.exp(-2*t) * Math.sin(y); |
117 |
dl |
1.1 |
} |
118 |
jsr166 |
1.2 |
static final double solu(double x, double y, double t) { |
119 |
|
|
return Math.exp(-2*t) * Math.sin(x) * Math.sin(y); |
120 |
dl |
1.1 |
} |
121 |
|
|
|
122 |
|
|
|
123 |
|
|
|
124 |
|
|
|
125 |
|
|
static final class Compute extends RecursiveAction { |
126 |
|
|
final int lb; |
127 |
|
|
final int ub; |
128 |
|
|
final int time; |
129 |
|
|
|
130 |
|
|
Compute(int lowerBound, int upperBound, int timestep) { |
131 |
|
|
lb = lowerBound; |
132 |
|
|
ub = upperBound; |
133 |
|
|
time = timestep; |
134 |
|
|
} |
135 |
jsr166 |
1.2 |
|
136 |
dl |
1.1 |
public void compute() { |
137 |
|
|
if (ub - lb > leafmaxcol) { |
138 |
|
|
int mid = (lb + ub) >>> 1; |
139 |
|
|
Compute left = new Compute(lb, mid, time); |
140 |
|
|
left.fork(); |
141 |
|
|
new Compute(mid, ub, time).compute(); |
142 |
|
|
left.join(); |
143 |
|
|
} |
144 |
|
|
else if (time == 0) // if first pass, initialize cells |
145 |
|
|
init(); |
146 |
|
|
else if (time %2 != 0) // alternate new/old |
147 |
|
|
compstripe(newm, oldm); |
148 |
|
|
else |
149 |
|
|
compstripe(oldm, newm); |
150 |
|
|
} |
151 |
|
|
|
152 |
|
|
|
153 |
|
|
/** Update all cells **/ |
154 |
|
|
final void compstripe(double[][] newMat, double[][] oldMat) { |
155 |
|
|
|
156 |
|
|
// manually mangled to reduce array indexing |
157 |
|
|
|
158 |
|
|
final int llb = (lb == 0) ? 1 : lb; |
159 |
|
|
final int lub = (ub == nx) ? nx - 1 : ub; |
160 |
|
|
|
161 |
|
|
double[] west; |
162 |
|
|
double[] row = oldMat[llb-1]; |
163 |
|
|
double[] east = oldMat[llb]; |
164 |
|
|
|
165 |
|
|
for (int a = llb; a < lub; a++) { |
166 |
|
|
|
167 |
|
|
west = row; |
168 |
|
|
row = east; |
169 |
|
|
east = oldMat[a+1]; |
170 |
|
|
|
171 |
|
|
double prev; |
172 |
|
|
double cell = row[0]; |
173 |
|
|
double next = row[1]; |
174 |
|
|
|
175 |
|
|
double[] nv = newMat[a]; |
176 |
|
|
|
177 |
|
|
for (int b = 1; b < ny-1; b++) { |
178 |
|
|
|
179 |
|
|
prev = cell; |
180 |
|
|
cell = next; |
181 |
|
|
double twoc = 2 * cell; |
182 |
|
|
next = row[b+1]; |
183 |
|
|
|
184 |
|
|
nv[b] = cell |
185 |
|
|
+ dtdysq * (prev - twoc + next) |
186 |
|
|
+ dtdxsq * (east[b] - twoc + west[b]); |
187 |
|
|
|
188 |
|
|
} |
189 |
|
|
} |
190 |
|
|
|
191 |
|
|
edges(newMat, llb, lub, tu + time * dt); |
192 |
|
|
} |
193 |
|
|
|
194 |
|
|
|
195 |
|
|
// the original version from cilk |
196 |
|
|
final void origcompstripe(double[][] newMat, double[][] oldMat) { |
197 |
jsr166 |
1.2 |
|
198 |
dl |
1.1 |
final int llb = (lb == 0) ? 1 : lb; |
199 |
|
|
final int lub = (ub == nx) ? nx - 1 : ub; |
200 |
|
|
|
201 |
|
|
for (int a = llb; a < lub; a++) { |
202 |
|
|
for (int b = 1; b < ny-1; b++) { |
203 |
|
|
double cell = oldMat[a][b]; |
204 |
|
|
double twoc = 2 * cell; |
205 |
|
|
newMat[a][b] = cell |
206 |
|
|
+ dtdxsq * (oldMat[a+1][b] - twoc + oldMat[a-1][b]) |
207 |
|
|
+ dtdysq * (oldMat[a][b+1] - twoc + oldMat[a][b-1]); |
208 |
|
|
|
209 |
|
|
} |
210 |
|
|
} |
211 |
|
|
|
212 |
|
|
edges(newMat, llb, lub, tu + time * dt); |
213 |
|
|
} |
214 |
|
|
|
215 |
|
|
|
216 |
|
|
/** Initialize all cells **/ |
217 |
|
|
final void init() { |
218 |
|
|
final int llb = (lb == 0) ? 1 : lb; |
219 |
|
|
final int lub = (ub == nx) ? nx - 1 : ub; |
220 |
|
|
|
221 |
|
|
for (int a = llb; a < lub; a++) { /* inner nodes */ |
222 |
|
|
double[] ov = oldm[a]; |
223 |
|
|
double x = xu + a * dx; |
224 |
|
|
double y = yu; |
225 |
|
|
for (int b = 1; b < ny-1; b++) { |
226 |
|
|
y += dy; |
227 |
|
|
ov[b] = f(x, y); |
228 |
|
|
} |
229 |
|
|
} |
230 |
|
|
|
231 |
|
|
edges(oldm, llb, lub, 0); |
232 |
|
|
|
233 |
|
|
} |
234 |
|
|
|
235 |
|
|
/** Fill in edges with boundary values **/ |
236 |
|
|
final void edges(double [][] m, int llb, int lub, double t) { |
237 |
|
|
|
238 |
|
|
for (int a = llb; a < lub; a++) { |
239 |
|
|
double[] v = m[a]; |
240 |
|
|
double x = xu + a * dx; |
241 |
|
|
v[0] = randa(x, t); |
242 |
|
|
v[ny-1] = randb(x, t); |
243 |
|
|
} |
244 |
|
|
|
245 |
|
|
if (lb == 0) { |
246 |
|
|
double[] v = m[0]; |
247 |
|
|
double y = yu; |
248 |
|
|
for (int b = 0; b < ny; b++) { |
249 |
|
|
y += dy; |
250 |
|
|
v[b] = randc(y, t); |
251 |
|
|
} |
252 |
|
|
} |
253 |
|
|
|
254 |
|
|
if (ub == nx) { |
255 |
jsr166 |
1.2 |
double[] v = m[nx - 1]; |
256 |
dl |
1.1 |
double y = yu; |
257 |
|
|
for (int b = 0; b < ny; b++) { |
258 |
|
|
y += dy; |
259 |
|
|
v[b] = randd(y, t); |
260 |
|
|
} |
261 |
|
|
} |
262 |
|
|
} |
263 |
|
|
} |
264 |
|
|
} |