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 |
jsr166 |
1.6 |
* http://creativecommons.org/publicdomain/zero/1.0/ |
5 |
dl |
1.1 |
*/ |
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.5 |
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 |
jsr166 |
1.9 |
/** the function being applied across the cells */ |
99 |
jsr166 |
1.2 |
static final double f(double x, double y) { |
100 |
|
|
return Math.sin(x) * Math.sin(y); |
101 |
dl |
1.1 |
} |
102 |
|
|
|
103 |
|
|
// random starting values |
104 |
|
|
|
105 |
jsr166 |
1.2 |
static final double randa(double x, double t) { |
106 |
|
|
return 0.0; |
107 |
dl |
1.1 |
} |
108 |
jsr166 |
1.2 |
static final double randb(double x, double t) { |
109 |
|
|
return Math.exp(-2*t) * Math.sin(x); |
110 |
dl |
1.1 |
} |
111 |
jsr166 |
1.2 |
static final double randc(double y, double t) { |
112 |
|
|
return 0.0; |
113 |
dl |
1.1 |
} |
114 |
jsr166 |
1.2 |
static final double randd(double y, double t) { |
115 |
|
|
return Math.exp(-2*t) * Math.sin(y); |
116 |
dl |
1.1 |
} |
117 |
jsr166 |
1.2 |
static final double solu(double x, double y, double t) { |
118 |
|
|
return Math.exp(-2*t) * Math.sin(x) * Math.sin(y); |
119 |
dl |
1.1 |
} |
120 |
|
|
|
121 |
|
|
static final class Compute extends RecursiveAction { |
122 |
|
|
final int lb; |
123 |
|
|
final int ub; |
124 |
|
|
final int time; |
125 |
|
|
|
126 |
|
|
Compute(int lowerBound, int upperBound, int timestep) { |
127 |
|
|
lb = lowerBound; |
128 |
|
|
ub = upperBound; |
129 |
|
|
time = timestep; |
130 |
|
|
} |
131 |
jsr166 |
1.2 |
|
132 |
dl |
1.1 |
public void compute() { |
133 |
|
|
if (ub - lb > leafmaxcol) { |
134 |
|
|
int mid = (lb + ub) >>> 1; |
135 |
|
|
Compute left = new Compute(lb, mid, time); |
136 |
|
|
left.fork(); |
137 |
|
|
new Compute(mid, ub, time).compute(); |
138 |
|
|
left.join(); |
139 |
|
|
} |
140 |
|
|
else if (time == 0) // if first pass, initialize cells |
141 |
|
|
init(); |
142 |
|
|
else if (time %2 != 0) // alternate new/old |
143 |
|
|
compstripe(newm, oldm); |
144 |
|
|
else |
145 |
|
|
compstripe(oldm, newm); |
146 |
|
|
} |
147 |
|
|
|
148 |
jsr166 |
1.4 |
/** Updates all cells. */ |
149 |
dl |
1.1 |
final void compstripe(double[][] newMat, double[][] oldMat) { |
150 |
|
|
|
151 |
|
|
// manually mangled to reduce array indexing |
152 |
|
|
|
153 |
|
|
final int llb = (lb == 0) ? 1 : lb; |
154 |
|
|
final int lub = (ub == nx) ? nx - 1 : ub; |
155 |
|
|
|
156 |
|
|
double[] west; |
157 |
|
|
double[] row = oldMat[llb-1]; |
158 |
|
|
double[] east = oldMat[llb]; |
159 |
|
|
|
160 |
|
|
for (int a = llb; a < lub; a++) { |
161 |
|
|
|
162 |
|
|
west = row; |
163 |
|
|
row = east; |
164 |
|
|
east = oldMat[a+1]; |
165 |
|
|
|
166 |
|
|
double prev; |
167 |
|
|
double cell = row[0]; |
168 |
|
|
double next = row[1]; |
169 |
|
|
|
170 |
|
|
double[] nv = newMat[a]; |
171 |
|
|
|
172 |
|
|
for (int b = 1; b < ny-1; b++) { |
173 |
|
|
|
174 |
|
|
prev = cell; |
175 |
|
|
cell = next; |
176 |
|
|
double twoc = 2 * cell; |
177 |
|
|
next = row[b+1]; |
178 |
|
|
|
179 |
|
|
nv[b] = cell |
180 |
|
|
+ dtdysq * (prev - twoc + next) |
181 |
|
|
+ dtdxsq * (east[b] - twoc + west[b]); |
182 |
|
|
} |
183 |
|
|
} |
184 |
|
|
|
185 |
jsr166 |
1.7 |
edges(newMat, llb, lub, tu + time * dt); |
186 |
dl |
1.1 |
} |
187 |
|
|
|
188 |
jsr166 |
1.9 |
/** the original version from cilk */ |
189 |
dl |
1.1 |
final void origcompstripe(double[][] newMat, double[][] oldMat) { |
190 |
jsr166 |
1.2 |
|
191 |
dl |
1.1 |
final int llb = (lb == 0) ? 1 : lb; |
192 |
|
|
final int lub = (ub == nx) ? nx - 1 : ub; |
193 |
|
|
|
194 |
|
|
for (int a = llb; a < lub; a++) { |
195 |
|
|
for (int b = 1; b < ny-1; b++) { |
196 |
|
|
double cell = oldMat[a][b]; |
197 |
|
|
double twoc = 2 * cell; |
198 |
|
|
newMat[a][b] = cell |
199 |
|
|
+ dtdxsq * (oldMat[a+1][b] - twoc + oldMat[a-1][b]) |
200 |
|
|
+ dtdysq * (oldMat[a][b+1] - twoc + oldMat[a][b-1]); |
201 |
|
|
} |
202 |
|
|
} |
203 |
|
|
|
204 |
|
|
edges(newMat, llb, lub, tu + time * dt); |
205 |
|
|
} |
206 |
|
|
|
207 |
jsr166 |
1.4 |
/** Initializes all cells. */ |
208 |
dl |
1.1 |
final void init() { |
209 |
|
|
final int llb = (lb == 0) ? 1 : lb; |
210 |
|
|
final int lub = (ub == nx) ? nx - 1 : ub; |
211 |
|
|
|
212 |
jsr166 |
1.3 |
for (int a = llb; a < lub; a++) { /* inner nodes */ |
213 |
dl |
1.1 |
double[] ov = oldm[a]; |
214 |
|
|
double x = xu + a * dx; |
215 |
|
|
double y = yu; |
216 |
|
|
for (int b = 1; b < ny-1; b++) { |
217 |
|
|
y += dy; |
218 |
|
|
ov[b] = f(x, y); |
219 |
|
|
} |
220 |
|
|
} |
221 |
|
|
|
222 |
|
|
edges(oldm, llb, lub, 0); |
223 |
|
|
} |
224 |
|
|
|
225 |
jsr166 |
1.4 |
/** Fills in edges with boundary values. */ |
226 |
dl |
1.1 |
final void edges(double [][] m, int llb, int lub, double t) { |
227 |
|
|
|
228 |
|
|
for (int a = llb; a < lub; a++) { |
229 |
|
|
double[] v = m[a]; |
230 |
|
|
double x = xu + a * dx; |
231 |
|
|
v[0] = randa(x, t); |
232 |
|
|
v[ny-1] = randb(x, t); |
233 |
|
|
} |
234 |
|
|
|
235 |
|
|
if (lb == 0) { |
236 |
|
|
double[] v = m[0]; |
237 |
|
|
double y = yu; |
238 |
|
|
for (int b = 0; b < ny; b++) { |
239 |
|
|
y += dy; |
240 |
|
|
v[b] = randc(y, t); |
241 |
|
|
} |
242 |
|
|
} |
243 |
|
|
|
244 |
|
|
if (ub == nx) { |
245 |
jsr166 |
1.2 |
double[] v = m[nx - 1]; |
246 |
dl |
1.1 |
double y = yu; |
247 |
|
|
for (int b = 0; b < ny; b++) { |
248 |
|
|
y += dy; |
249 |
|
|
v[b] = randd(y, t); |
250 |
|
|
} |
251 |
|
|
} |
252 |
|
|
} |
253 |
|
|
} |
254 |
|
|
} |