test/loops/Heat.java

/*
 * Written by Doug Lea with assistance from members of JCP JSR-166
 * Expert Group and released to the public domain, as explained at
 * http://creativecommons.org/publicdomain/zero/1.0/
 */

// Adapted from a cilk benchmark

import java.util.concurrent.*;

public class Heat {
    static final long NPS = (1000L * 1000 * 1000);

    // Parameters
    static int nx;
    static int ny;
    static int nt;
    static int leafmaxcol;

    // the matrix representing the cells
    static double[][] newm;

    // alternating workspace matrix
    static double[][] oldm;

    public static void main(String[] args) throws Exception {
        int procs = 0;
        nx = 4096;
        ny = 1024;
        nt = 1000;
        leafmaxcol = 16;

        try {
            if (args.length > 0)
                procs = Integer.parseInt(args[0]);
            if (args.length > 1)
                nx = Integer.parseInt(args[1]);
            if (args.length > 2)
                ny = Integer.parseInt(args[2]);
            if (args.length > 3)
                nt = Integer.parseInt(args[3]);
            if (args.length > 4)
                leafmaxcol = Integer.parseInt(args[4]);
        }
        catch (Exception e) {
            System.out.println("Usage: java Heat threads rows cols steps granularity");
            return;
        }

        ForkJoinPool g = (procs == 0) ? new ForkJoinPool() :
            new ForkJoinPool(procs);

        System.out.print("parallelism = " + g.getParallelism());
        System.out.print(" granularity = " + leafmaxcol);
        System.out.print(" rows = " + nx);
        System.out.print(" columns = " + ny);
        System.out.println(" steps = " + nt);

        oldm = new double[nx][ny];
        newm = new double[nx][ny];

        for (int i = 0; i < 5; ++i) {
            long last = System.nanoTime();
            RecursiveAction main = new RecursiveAction() {
                    public void compute() {
                        for (int timestep = 0; timestep <= nt; timestep++) {
                            (new Compute(0, nx, timestep)).invoke();
                        }
                    }
                };
            g.invoke(main);
            double elapsed = elapsedTime(last);
            System.out.printf("time: %7.3f", elapsed);
            System.out.println();
        }
        System.out.println(g);
        g.shutdown();
    }

    static double elapsedTime(long startTime) {
        return (double)(System.nanoTime() - startTime) / NPS;
    }

    // constants (at least for this demo)
    static final double xu = 0.0;
    static final double xo = 1.570796326794896558;
    static final double yu = 0.0;
    static final double yo = 1.570796326794896558;
    static final double tu = 0.0;
    static final double to = 0.0000001;

    static final double dx = (xo - xu) / (nx - 1);
    static final double dy = (yo - yu) / (ny - 1);
    static final double dt = (to - tu) / nt;
    static final double dtdxsq = dt / (dx * dx);
    static final double dtdysq = dt / (dy * dy);


    // the function being applied across the cells
    static final double f(double x, double y) {
        return Math.sin(x) * Math.sin(y);
    }

    // random starting values

    static final double randa(double x, double t) {
        return 0.0;
    }
    static final double randb(double x, double t) {
        return Math.exp(-2*t) * Math.sin(x);
    }
    static final double randc(double y, double t) {
        return 0.0;
    }
    static final double randd(double y, double t) {
        return Math.exp(-2*t) * Math.sin(y);
    }
    static final double solu(double x, double y, double t) {
        return Math.exp(-2*t) * Math.sin(x) * Math.sin(y);
    }


    static final class Compute extends RecursiveAction {
        final int lb;
        final int ub;
        final int time;

        Compute(int lowerBound, int upperBound, int timestep) {
            lb = lowerBound;
            ub = upperBound;
            time = timestep;
        }

        public void compute() {
            if (ub - lb > leafmaxcol) {
                int mid = (lb + ub) >>> 1;
                Compute left = new Compute(lb, mid, time);
                left.fork();
                new Compute(mid, ub, time).compute();
                left.join();
            }
            else if (time == 0)     // if first pass, initialize cells
                init();
            else if (time %2 != 0)  // alternate new/old
                compstripe(newm, oldm);
            else
                compstripe(oldm, newm);
        }


        /** Updates all cells. */
        final void compstripe(double[][] newMat, double[][] oldMat) {

            // manually mangled to reduce array indexing

            final int llb = (lb == 0)  ? 1 : lb;
            final int lub = (ub == nx) ? nx - 1 : ub;

            double[] west;
            double[] row = oldMat[llb-1];
            double[] east = oldMat[llb];

            for (int a = llb; a < lub; a++) {

                west = row;
                row =  east;
                east = oldMat[a+1];

                double prev;
                double cell = row[0];
                double next = row[1];

                double[] nv = newMat[a];

                for (int b = 1; b < ny-1; b++) {

                    prev = cell;
                    cell = next;
                    double twoc = 2 * cell;
                    next = row[b+1];

                    nv[b] = cell
                        + dtdysq * (prev    - twoc + next)
                        + dtdxsq * (east[b] - twoc + west[b]);

                }
            }

            edges(newMat, llb, lub,  tu + time * dt);
        }


        // the original version from cilk
        final void origcompstripe(double[][] newMat, double[][] oldMat) {

            final int llb = (lb == 0)  ? 1 : lb;
            final int lub = (ub == nx) ? nx - 1 : ub;

            for (int a = llb; a < lub; a++) {
                for (int b = 1; b < ny-1; b++) {
                    double cell = oldMat[a][b];
                    double twoc = 2 * cell;
                    newMat[a][b] = cell
                        + dtdxsq * (oldMat[a+1][b] - twoc + oldMat[a-1][b])
                        + dtdysq * (oldMat[a][b+1] - twoc + oldMat[a][b-1]);

                }
            }

            edges(newMat, llb, lub,  tu + time * dt);
        }


        /** Initializes all cells. */
        final void init() {
            final int llb = (lb == 0) ? 1 : lb;
            final int lub = (ub == nx) ? nx - 1 : ub;

            for (int a = llb; a < lub; a++) {   /* inner nodes */
                double[] ov = oldm[a];
                double x = xu + a * dx;
                double y = yu;
                for (int b = 1; b < ny-1; b++) {
                    y += dy;
                    ov[b] = f(x, y);
                }
            }

            edges(oldm, llb, lub, 0);

        }

        /** Fills in edges with boundary values. */
        final void edges(double [][] m, int llb, int lub, double t) {

            for (int a = llb; a < lub; a++) {
                double[] v = m[a];
                double x = xu + a * dx;
                v[0] = randa(x, t);
                v[ny-1] = randb(x, t);
            }

            if (lb == 0) {
                double[] v = m[0];
                double y = yu;
                for (int b = 0; b < ny; b++) {
                    y += dy;
                    v[b] = randc(y, t);
                }
            }

            if (ub == nx) {
                double[] v = m[nx - 1];
                double y = yu;
                for (int b = 0; b < ny; b++) {
                    y += dy;
                    v[b] = randd(y, t);
                }
            }
        }
    }
}
Revision:	1.6
Committed:	Tue Mar 15 19:47:05 2011 UTC (13 years, 2 months ago) by jsr166
Branch:	MAIN
CVS Tags:	release-1_7_0
Changes since 1.5:	+1 -1 lines
Log Message:	Update Creative Commons license URL in legal notices
#	Content
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(-2t) 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(-2t) Math.sin(y);
117	}
118	static final double solu(double x, double y, double t) {
119	return Math.exp(-2t) 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
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
198	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	/** Initializes 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	/** Fills 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	double[] v = m[nx - 1];
256	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	}