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, 8 months ago) by jsr166
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.8:	+2 -9 lines
Log Message:	style nits
#	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	jsr166	1.6	* http://creativecommons.org/publicdomain/zero/1.0/
5	dl	1.1	*/
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.5	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	jsr166	1.9	/** the function being applied across the cells */
99	jsr166	1.2	static final double f(double x, double y) {
100			return Math.sin(x) * Math.sin(y);
101	dl	1.1	}
102
103			// random starting values
104
105	jsr166	1.2	static final double randa(double x, double t) {
106			return 0.0;
107	dl	1.1	}
108	jsr166	1.2	static final double randb(double x, double t) {
109			return Math.exp(-2t) Math.sin(x);
110	dl	1.1	}
111	jsr166	1.2	static final double randc(double y, double t) {
112			return 0.0;
113	dl	1.1	}
114	jsr166	1.2	static final double randd(double y, double t) {
115			return Math.exp(-2t) Math.sin(y);
116	dl	1.1	}
117	jsr166	1.2	static final double solu(double x, double y, double t) {
118			return Math.exp(-2t) Math.sin(x) * Math.sin(y);
119	dl	1.1	}
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	jsr166	1.2
132	dl	1.1	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	jsr166	1.4	/** Updates all cells. */
149	dl	1.1	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	jsr166	1.7	edges(newMat, llb, lub, tu + time * dt);
186	dl	1.1	}
187
188	jsr166	1.9	/** the original version from cilk */
189	dl	1.1	final void origcompstripe(double[][] newMat, double[][] oldMat) {
190	jsr166	1.2
191	dl	1.1	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	jsr166	1.4	/** Initializes all cells. */
208	dl	1.1	final void init() {
209			final int llb = (lb == 0) ? 1 : lb;
210			final int lub = (ub == nx) ? nx - 1 : ub;
211
212	jsr166	1.3	for (int a = llb; a < lub; a++) { /* inner nodes */
213	dl	1.1	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	jsr166	1.4	/** Fills in edges with boundary values. */
226	dl	1.1	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	jsr166	1.2	double[] v = m[nx - 1];
246	dl	1.1	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			}