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 |
/** the function being applied across the cells */ |
99 |
static final double f(double x, double y) { |
100 |
return Math.sin(x) * Math.sin(y); |
101 |
} |
102 |
|
103 |
// random starting values |
104 |
|
105 |
static final double randa(double x, double t) { |
106 |
return 0.0; |
107 |
} |
108 |
static final double randb(double x, double t) { |
109 |
return Math.exp(-2*t) * Math.sin(x); |
110 |
} |
111 |
static final double randc(double y, double t) { |
112 |
return 0.0; |
113 |
} |
114 |
static final double randd(double y, double t) { |
115 |
return Math.exp(-2*t) * Math.sin(y); |
116 |
} |
117 |
static final double solu(double x, double y, double t) { |
118 |
return Math.exp(-2*t) * Math.sin(x) * Math.sin(y); |
119 |
} |
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 |
|
132 |
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 |
/** Updates all cells. */ |
149 |
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 |
edges(newMat, llb, lub, tu + time * dt); |
186 |
} |
187 |
|
188 |
/** the original version from cilk */ |
189 |
final void origcompstripe(double[][] newMat, double[][] oldMat) { |
190 |
|
191 |
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 |
/** Initializes all cells. */ |
208 |
final void init() { |
209 |
final int llb = (lb == 0) ? 1 : lb; |
210 |
final int lub = (ub == nx) ? nx - 1 : ub; |
211 |
|
212 |
for (int a = llb; a < lub; a++) { /* inner nodes */ |
213 |
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 |
/** Fills in edges with boundary values. */ |
226 |
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 |
double[] v = m[nx - 1]; |
246 |
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 |
} |