test/loops/IntegrateGamma.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
 */

import java.util.concurrent.*;

/**
 * Adapted from FJTask version.
 * Sample program using Guassian Quadrature for numerical integration.
 * Inspired by a
 * <A href="http://www.cs.uga.edu/~dkl/filaments/dist.html"> Filaments</A>
 * demo program.
 *
 */

public class IntegrateGamma {
    /** for time conversion */
    static final long NPS = (1000L * 1000 * 1000);
    public static void main(String[] args) {
        int procs = 0;
        double start = 1.0;
        double end = 96.0;
        int exp = 5;
        try {
            if (args.length > 0)
                procs = Integer.parseInt(args[0]);
            if (args.length > 1)
                start = new Double(args[1]).doubleValue();
            if (args.length > 2)
                end = new Double(args[2]).doubleValue();
            if (args.length > 3)
                exp = Integer.parseInt(args[3]);
        }
        catch (Exception e) {
            System.out.println("Usage: java IntegrateGamma <threads> <lower bound> <upper bound> <exponent>\n (for example 2 1 48 5).");
            return;
        }

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

        System.out.println("Integrating from " + start + " to " + end + " exponent: " + exp + " parallelism " + g.getParallelism());

        Function f = new SampleFunction(exp);
        for (int i = 0; i < 10; ++i) {
            Integrator integrator = new Integrator(f, 0.001, g);
            long last = System.nanoTime();
            double result = integrator.integral(start, end);
            double elapsed = elapsedTime(last);
            System.out.printf("time: %7.3f", elapsed);
            System.out.println(" Answer = " + result);
        }
        System.out.println(g);
        g.shutdown();
    }

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

    /*
      This is all set up as if it were part of a more serious
      framework, but is for now just a demo, with all
      classes declared as static within Integrate
    */

    /** A function to be integrated */
    static interface Function {
        double compute(double x);
    }

    /**
     * Sample from filaments demo.
     * Computes (2*n-1)*(x^(2*n-1)) for all odd values.
     */
    static class SampleFunction implements Function {
        final int n;
        SampleFunction(int n) { this.n = n; }

        public double compute(double x) {
            double power = x;
            double xsq = x * x;
            double val = power;
            double di = 1.0;
            for (int i = n - 1; i > 0; --i) {
                di += 2.0;
                power *= xsq;
                val += di * power;
            }
            return val;
        }
    }


    static class Integrator {
        final Function f;      // The function to integrate
        final double errorTolerance;
        final ForkJoinPool g;

        Integrator(Function f, double errorTolerance, ForkJoinPool g) {
            this.f = f;
            this.errorTolerance = errorTolerance;
            this.g = g;
        }

        double integral(double lowerBound, double upperBound) {
            double f_lower = f.compute(lowerBound);
            double f_upper = f.compute(upperBound);
            double initialArea = 0.5 * (upperBound-lowerBound) * (f_upper + f_lower);
            Quad q = new Quad(lowerBound, upperBound,
                              f_lower, f_upper,
                              initialArea);
            g.invoke(q);
            return q.area;
        }


        /**
         * FJTask to recursively perform the quadrature.
         * Algorithm:
         *  Compute the area from lower bound to the center point of interval,
         *  and from the center point to the upper bound. If this
         *  differs from the value from lower to upper by more than
         *  the error tolerance, recurse on each half.
         */
        final class Quad extends RecursiveAction {
            final double left;       // lower bound
            final double right;      // upper bound
            final double f_left;     // value of the function evaluated at left
            final double f_right;    // value of the function evaluated at right

            // Area initialized with original estimate from left to right.
            // It is replaced with refined value.
            volatile double area;

            Quad(double left, double right,
                 double f_left, double f_right,
                 double area) {
                this.left = left;
                this.right = right;
                this.f_left = f_left;
                this.f_right = f_right;
                this.area = area;
            }

            public void compute() {
                double center = 0.5 * (left + right);
                double f_center = f.compute(center);

                double leftArea  = 0.5 * (center - left)  * (f_left + f_center);
                double rightArea = 0.5 * (right - center) * (f_center + f_right);
                double sum = leftArea + rightArea;

                double diff = sum - area;
                if (diff < 0) diff = -diff;

                if (diff >= errorTolerance) {
                    Quad q1 = new Quad(left,   center, f_left,   f_center, leftArea);
                    q1.fork();
                    Quad q2 = new Quad(center, right,  f_center, f_right,  rightArea);
                    q2.compute();
                    q1.join();
                    sum = q1.area + q2.area;
                }

                area = sum;
            }
        }
    }

}


Revision:	1.3
Committed:	Sat Oct 16 16:22:57 2010 UTC (13 years, 6 months ago) by jsr166
Branch:	MAIN
Changes since 1.2:	+5 -5 lines
Log Message:	coding style
#	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			import java.util.concurrent.*;
8
9			/**
10			* Adapted from FJTask version.
11			* Sample program using Guassian Quadrature for numerical integration.
12	jsr166	1.2	* Inspired by a
13	dl	1.1	* <A href="http://www.cs.uga.edu/~dkl/filaments/dist.html"> Filaments</A>
14			* demo program.
15	jsr166	1.2	*
16	dl	1.1	*/
17
18			public class IntegrateGamma {
19			/** for time conversion */
20			static final long NPS = (1000L * 1000 * 1000);
21			public static void main(String[] args) {
22			int procs = 0;
23			double start = 1.0;
24			double end = 96.0;
25			int exp = 5;
26			try {
27			if (args.length > 0)
28			procs = Integer.parseInt(args[0]);
29			if (args.length > 1)
30			start = new Double(args[1]).doubleValue();
31			if (args.length > 2)
32			end = new Double(args[2]).doubleValue();
33	jsr166	1.2	if (args.length > 3)
34	dl	1.1	exp = Integer.parseInt(args[3]);
35			}
36			catch (Exception e) {
37			System.out.println("Usage: java IntegrateGamma <threads> <lower bound> <upper bound> <exponent>\n (for example 2 1 48 5).");
38			return;
39			}
40
41	jsr166	1.2	ForkJoinPool g = procs == 0? new ForkJoinPool() :
42	dl	1.1	new ForkJoinPool(procs);
43
44			System.out.println("Integrating from " + start + " to " + end + " exponent: " + exp + " parallelism " + g.getParallelism());
45	jsr166	1.2
46	dl	1.1	Function f = new SampleFunction(exp);
47			for (int i = 0; i < 10; ++i) {
48			Integrator integrator = new Integrator(f, 0.001, g);
49			long last = System.nanoTime();
50			double result = integrator.integral(start, end);
51			double elapsed = elapsedTime(last);
52			System.out.printf("time: %7.3f", elapsed);
53			System.out.println(" Answer = " + result);
54	jsr166	1.2	}
55	dl	1.1	System.out.println(g);
56			g.shutdown();
57			}
58
59			static double elapsedTime(long startTime) {
60			return (double)(System.nanoTime() - startTime) / NPS;
61			}
62
63			/*
64			This is all set up as if it were part of a more serious
65			framework, but is for now just a demo, with all
66			classes declared as static within Integrate
67			*/
68
69	jsr166	1.3	/** A function to be integrated */
70	dl	1.1	static interface Function {
71			double compute(double x);
72			}
73
74			/**
75			* Sample from filaments demo.
76	jsr166	1.3	* Computes (2n-1)(x^(2*n-1)) for all odd values.
77			*/
78	dl	1.1	static class SampleFunction implements Function {
79			final int n;
80			SampleFunction(int n) { this.n = n; }
81
82	jsr166	1.3	public double compute(double x) {
83	dl	1.1	double power = x;
84			double xsq = x * x;
85			double val = power;
86			double di = 1.0;
87			for (int i = n - 1; i > 0; --i) {
88			di += 2.0;
89			power *= xsq;
90			val += di * power;
91			}
92			return val;
93			}
94			}
95
96
97			static class Integrator {
98			final Function f; // The function to integrate
99			final double errorTolerance;
100			final ForkJoinPool g;
101
102			Integrator(Function f, double errorTolerance, ForkJoinPool g) {
103			this.f = f;
104			this.errorTolerance = errorTolerance;
105			this.g = g;
106			}
107
108			double integral(double lowerBound, double upperBound) {
109			double f_lower = f.compute(lowerBound);
110			double f_upper = f.compute(upperBound);
111			double initialArea = 0.5 * (upperBound-lowerBound) * (f_upper + f_lower);
112			Quad q = new Quad(lowerBound, upperBound,
113			f_lower, f_upper,
114			initialArea);
115			g.invoke(q);
116			return q.area;
117			}
118
119	jsr166	1.2
120			/**
121	dl	1.1	* FJTask to recursively perform the quadrature.
122			* Algorithm:
123			* Compute the area from lower bound to the center point of interval,
124			* and from the center point to the upper bound. If this
125			* differs from the value from lower to upper by more than
126			* the error tolerance, recurse on each half.
127	jsr166	1.3	*/
128	dl	1.1	final class Quad extends RecursiveAction {
129			final double left; // lower bound
130			final double right; // upper bound
131			final double f_left; // value of the function evaluated at left
132			final double f_right; // value of the function evaluated at right
133	jsr166	1.2
134	dl	1.1	// Area initialized with original estimate from left to right.
135			// It is replaced with refined value.
136			volatile double area;
137	jsr166	1.2
138			Quad(double left, double right,
139			double f_left, double f_right,
140	dl	1.1	double area) {
141			this.left = left;
142			this.right = right;
143			this.f_left = f_left;
144			this.f_right = f_right;
145			this.area = area;
146			}
147	jsr166	1.2
148	dl	1.1	public void compute() {
149			double center = 0.5 * (left + right);
150	jsr166	1.2	double f_center = f.compute(center);
151
152			double leftArea = 0.5 * (center - left) * (f_left + f_center);
153	dl	1.1	double rightArea = 0.5 * (right - center) * (f_center + f_right);
154			double sum = leftArea + rightArea;
155	jsr166	1.2
156	dl	1.1	double diff = sum - area;
157			if (diff < 0) diff = -diff;
158	jsr166	1.2
159			if (diff >= errorTolerance) {
160	dl	1.1	Quad q1 = new Quad(left, center, f_left, f_center, leftArea);
161			q1.fork();
162			Quad q2 = new Quad(center, right, f_center, f_right, rightArea);
163			q2.compute();
164			q1.join();
165	jsr166	1.2	sum = q1.area + q2.area;
166	dl	1.1	}
167	jsr166	1.2
168	dl	1.1	area = sum;
169			}
170			}
171			}
172
173			}
174
175	jsr166	1.2