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/publicdomain/zero/1.0/ |
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 |
ForkJoinPool g = (procs == 0) ? new ForkJoinPool() : |
51 |
new ForkJoinPool(procs); |
52 |
|
53 |
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 |
|
59 |
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 |
static final double dt = (to - tu) / nt; |
95 |
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 |
static final double f(double x, double y) { |
101 |
return Math.sin(x) * Math.sin(y); |
102 |
} |
103 |
|
104 |
// random starting values |
105 |
|
106 |
static final double randa(double x, double t) { |
107 |
return 0.0; |
108 |
} |
109 |
static final double randb(double x, double t) { |
110 |
return Math.exp(-2*t) * Math.sin(x); |
111 |
} |
112 |
static final double randc(double y, double t) { |
113 |
return 0.0; |
114 |
} |
115 |
static final double randd(double y, double t) { |
116 |
return Math.exp(-2*t) * Math.sin(y); |
117 |
} |
118 |
static final double solu(double x, double y, double t) { |
119 |
return Math.exp(-2*t) * Math.sin(x) * Math.sin(y); |
120 |
} |
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 |
|
136 |
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 |
/** Updates 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 |
edges(newMat, llb, lub, tu + time * dt); |
191 |
} |
192 |
|
193 |
|
194 |
// the original version from cilk |
195 |
final void origcompstripe(double[][] newMat, double[][] oldMat) { |
196 |
|
197 |
final int llb = (lb == 0) ? 1 : lb; |
198 |
final int lub = (ub == nx) ? nx - 1 : ub; |
199 |
|
200 |
for (int a = llb; a < lub; a++) { |
201 |
for (int b = 1; b < ny-1; b++) { |
202 |
double cell = oldMat[a][b]; |
203 |
double twoc = 2 * cell; |
204 |
newMat[a][b] = cell |
205 |
+ dtdxsq * (oldMat[a+1][b] - twoc + oldMat[a-1][b]) |
206 |
+ dtdysq * (oldMat[a][b+1] - twoc + oldMat[a][b-1]); |
207 |
} |
208 |
} |
209 |
|
210 |
edges(newMat, llb, lub, tu + time * dt); |
211 |
} |
212 |
|
213 |
|
214 |
/** Initializes all cells. */ |
215 |
final void init() { |
216 |
final int llb = (lb == 0) ? 1 : lb; |
217 |
final int lub = (ub == nx) ? nx - 1 : ub; |
218 |
|
219 |
for (int a = llb; a < lub; a++) { /* inner nodes */ |
220 |
double[] ov = oldm[a]; |
221 |
double x = xu + a * dx; |
222 |
double y = yu; |
223 |
for (int b = 1; b < ny-1; b++) { |
224 |
y += dy; |
225 |
ov[b] = f(x, y); |
226 |
} |
227 |
} |
228 |
|
229 |
edges(oldm, llb, lub, 0); |
230 |
} |
231 |
|
232 |
/** Fills in edges with boundary values. */ |
233 |
final void edges(double [][] m, int llb, int lub, double t) { |
234 |
|
235 |
for (int a = llb; a < lub; a++) { |
236 |
double[] v = m[a]; |
237 |
double x = xu + a * dx; |
238 |
v[0] = randa(x, t); |
239 |
v[ny-1] = randb(x, t); |
240 |
} |
241 |
|
242 |
if (lb == 0) { |
243 |
double[] v = m[0]; |
244 |
double y = yu; |
245 |
for (int b = 0; b < ny; b++) { |
246 |
y += dy; |
247 |
v[b] = randc(y, t); |
248 |
} |
249 |
} |
250 |
|
251 |
if (ub == nx) { |
252 |
double[] v = m[nx - 1]; |
253 |
double y = yu; |
254 |
for (int b = 0; b < ny; b++) { |
255 |
y += dy; |
256 |
v[b] = randd(y, t); |
257 |
} |
258 |
} |
259 |
} |
260 |
} |
261 |
} |