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/licenses/publicdomain
 */

// 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);
        }


        /** Update 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);
        }


        /** Initialize 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);

        }

        /** Fill 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.2
Committed:	Mon Sep 20 20:42:37 2010 UTC (13 years, 7 months ago) by jsr166
Branch:	MAIN
Changes since 1.1:	+19 -19 lines
Log Message:	whitespace
#	User	Rev	Content
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			* http://creativecommons.org/licenses/publicdomain
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	jsr166	1.2	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
99			// the function being applied across the cells
100	jsr166	1.2	static final double f(double x, double y) {
101			return Math.sin(x) * Math.sin(y);
102	dl	1.1	}
103
104			// random starting values
105
106	jsr166	1.2	static final double randa(double x, double t) {
107			return 0.0;
108	dl	1.1	}
109	jsr166	1.2	static final double randb(double x, double t) {
110			return Math.exp(-2t) Math.sin(x);
111	dl	1.1	}
112	jsr166	1.2	static final double randc(double y, double t) {
113			return 0.0;
114	dl	1.1	}
115	jsr166	1.2	static final double randd(double y, double t) {
116			return Math.exp(-2t) Math.sin(y);
117	dl	1.1	}
118	jsr166	1.2	static final double solu(double x, double y, double t) {
119			return Math.exp(-2t) Math.sin(x) * Math.sin(y);
120	dl	1.1	}
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	jsr166	1.2
136	dl	1.1	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			/ Update 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	jsr166	1.2
198	dl	1.1	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			/ Initialize 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			/ Fill 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	jsr166	1.2	double[] v = m[nx - 1];
256	dl	1.1	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			}