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.9
Committed:	Sat Sep 12 19:09:00 2015 UTC (8 years, 7 months ago) by jsr166
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.8:	+2 -9 lines
Log Message:	style nits
#	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	/** 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(-2t) 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(-2t) Math.sin(y);
116	}
117	static final double solu(double x, double y, double t) {
118	return Math.exp(-2t) 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	}