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

//import jsr166y.*;
import java.util.concurrent.*;

/**
 * Divide and Conquer matrix multiply demo
 */
public class MatrixMultiply {

    /** for time conversion */
    static final long NPS = (1000L * 1000 * 1000);

    static final int DEFAULT_GRANULARITY = 16; // 32;

    /**
     * The quadrant size at which to stop recursing down
     * and instead directly multiply the matrices.
     * Must be a power of two. Minimum value is 2.
     */
    static int granularity = DEFAULT_GRANULARITY;

    public static void main(String[] args) throws Exception {

        final String usage = "Usage: java MatrixMultiply <threads> <matrix size (must be a power of two)> [<granularity>] \n Size and granularity must be powers of two.\n For example, try java MatrixMultiply 2 512 16";

        int procs = 0;
        int n = 2048;
        int runs = 32;
        try {
            if (args.length > 0)
                procs = Integer.parseInt(args[0]);
            if (args.length > 1)
                n = Integer.parseInt(args[1]);
            if (args.length > 2)
                granularity = Integer.parseInt(args[2]);
            if (args.length > 3)
                runs = Integer.parseInt(args[2]);
        }

        catch (Exception e) {
            System.out.println(usage);
            return;
        }

        if ( ((n & (n - 1)) != 0) ||
             ((granularity & (granularity - 1)) != 0) ||
             granularity < 2) {
            System.out.println(usage);
            return;
        }

        ForkJoinPool pool = (procs == 0) ? ForkJoinPool.commonPool() :
            new ForkJoinPool(procs);
        System.out.println("procs: " + pool.getParallelism() +
                           " n: " + n + " granularity: " + granularity +
                           " runs: " + runs);

        float[][] a = new float[n][n];
        float[][] b = new float[n][n];
        float[][] c = new float[n][n];

        for (int i = 0; i < runs; ++i) {
            init(a, b, n);
            long start = System.nanoTime();
            new Multiplier(a, 0, 0, b, 0, 0, c, 0, 0, n).invoke();
            long time = System.nanoTime() - start;
            double secs = ((double)time) / NPS;
            Thread.sleep(100);
            System.out.printf("Time: %7.3f ", secs);
            if ((i & 3) == 3) System.out.println();
            // check(c, n);
        }
        System.out.println(pool.toString());
        if (pool != ForkJoinPool.commonPool())
            pool.shutdown();
        Thread.sleep(100);
    }

    // To simplify checking, fill with all 1's. Answer should be all n's.
    static void init(float[][] a, float[][] b, int n) {
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                a[i][j] = 1.0F;
                b[i][j] = 1.0F;
            }
        }
    }

    static void check(float[][] c, int n) {
        for (int i = 0; i < n; i++ ) {
            for (int j = 0; j < n; j++ ) {
                if (c[i][j] != n) {
                    throw new Error("Check Failed at [" + i +"]["+j+"]: " + c[i][j]);
                }
            }
        }
    }

    /**
     * Multiply matrices AxB by dividing into quadrants, using algorithm:
     * <pre>
     *      A      x      B
     *
     *  A11 | A12     B11 | B12     A11*B11 | A11*B12     A12*B21 | A12*B22
     * |----+----| x |----+----| = |--------+--------| + |---------+-------|
     *  A21 | A22     B21 | B21     A21*B11 | A21*B21     A22*B21 | A22*B22
     * </pre>
     */
    static class Multiplier extends RecursiveAction {
        final float[][] A;   // Matrix A
        final int aRow;      // first row    of current quadrant of A
        final int aCol;      // first column of current quadrant of A

        final float[][] B;   // Similarly for B
        final int bRow;
        final int bCol;

        final float[][] C;   // Similarly for result matrix C
        final int cRow;
        final int cCol;

        final int size;      // number of elements in current quadrant

        Multiplier(float[][] A, int aRow, int aCol,
                   float[][] B, int bRow, int bCol,
                   float[][] C, int cRow, int cCol,
                   int size) {
            this.A = A; this.aRow = aRow; this.aCol = aCol;
            this.B = B; this.bRow = bRow; this.bCol = bCol;
            this.C = C; this.cRow = cRow; this.cCol = cCol;
            this.size = size;
        }

        public void compute() {

            if (size <= granularity) {
                multiplyStride2();
            }

            else {
                int h = size / 2;

                invokeAll(new Seq2[] {
                    seq(new Multiplier(A, aRow,   aCol,    // A11
                                       B, bRow,   bCol,    // B11
                                       C, cRow,   cCol,    // C11
                                       h),
                        new Multiplier(A, aRow,   aCol+h,  // A12
                                       B, bRow+h, bCol,    // B21
                                       C, cRow,   cCol,    // C11
                                       h)),

                    seq(new Multiplier(A, aRow,   aCol,    // A11
                                       B, bRow,   bCol+h,  // B12
                                       C, cRow,   cCol+h,  // C12
                                       h),
                        new Multiplier(A, aRow,   aCol+h,  // A12
                                       B, bRow+h, bCol+h,  // B22
                                       C, cRow,   cCol+h,  // C12
                                       h)),

                    seq(new Multiplier(A, aRow+h, aCol,    // A21
                                       B, bRow,   bCol,    // B11
                                       C, cRow+h, cCol,    // C21
                                       h),
                        new Multiplier(A, aRow+h, aCol+h,  // A22
                                       B, bRow+h, bCol,    // B21
                                       C, cRow+h, cCol,    // C21
                                       h)),

                    seq(new Multiplier(A, aRow+h, aCol,    // A21
                                       B, bRow,   bCol+h,  // B12
                                       C, cRow+h, cCol+h,  // C22
                                       h),
                        new Multiplier(A, aRow+h, aCol+h,  // A22
                                       B, bRow+h, bCol+h,  // B22
                                       C, cRow+h, cCol+h,  // C22
                                       h))
                });
            }
        }

        /**
         * Version of matrix multiplication that steps 2 rows and columns
         * at a time. Adapted from Cilk demos.
         * Note that the results are added into C, not just set into C.
         * This works well here because Java array elements
         * are created with all zero values.
         */
        void multiplyStride2() {
            for (int j = 0; j < size; j+=2) {
                for (int i = 0; i < size; i +=2) {

                    float[] a0 = A[aRow+i];
                    float[] a1 = A[aRow+i+1];

                    float s00 = 0.0F;
                    float s01 = 0.0F;
                    float s10 = 0.0F;
                    float s11 = 0.0F;

                    for (int k = 0; k < size; k+=2) {

                        float[] b0 = B[bRow+k];

                        s00 += a0[aCol+k]   * b0[bCol+j];
                        s10 += a1[aCol+k]   * b0[bCol+j];
                        s01 += a0[aCol+k]   * b0[bCol+j+1];
                        s11 += a1[aCol+k]   * b0[bCol+j+1];

                        float[] b1 = B[bRow+k+1];

                        s00 += a0[aCol+k+1] * b1[bCol+j];
                        s10 += a1[aCol+k+1] * b1[bCol+j];
                        s01 += a0[aCol+k+1] * b1[bCol+j+1];
                        s11 += a1[aCol+k+1] * b1[bCol+j+1];
                    }

                    C[cRow+i]  [cCol+j]   += s00;
                    C[cRow+i]  [cCol+j+1] += s01;
                    C[cRow+i+1][cCol+j]   += s10;
                    C[cRow+i+1][cCol+j+1] += s11;
                }
            }
        }

    }

    static Seq2 seq(RecursiveAction task1,
                    RecursiveAction task2) {
        return new Seq2(task1, task2);
    }

    static final class Seq2 extends RecursiveAction {
        final RecursiveAction fst;
        final RecursiveAction snd;
        public Seq2(RecursiveAction task1, RecursiveAction task2) {
            fst = task1;
            snd = task2;
        }
        public void compute() {
            fst.invoke();
            snd.invoke();
        }
    }
}
Revision:	1.9
Committed:	Sat Sep 12 19:39:26 2015 UTC (8 years, 8 months ago) by dl
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.8:	+9 -6 lines
Log Message:	Use commonPool
#	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.5	* http://creativecommons.org/publicdomain/zero/1.0/
5	dl	1.1	*/
6
7			//import jsr166y.*;
8			import java.util.concurrent.*;
9
10			/**
11			* Divide and Conquer matrix multiply demo
12	jsr166	1.3	*/
13	dl	1.1	public class MatrixMultiply {
14
15			/** for time conversion */
16			static final long NPS = (1000L * 1000 * 1000);
17
18	dl	1.9	static final int DEFAULT_GRANULARITY = 16; // 32;
19	dl	1.1
20	jsr166	1.3	/**
21			* The quadrant size at which to stop recursing down
22	dl	1.1	* and instead directly multiply the matrices.
23			* Must be a power of two. Minimum value is 2.
24	jsr166	1.3	*/
25	dl	1.1	static int granularity = DEFAULT_GRANULARITY;
26
27			public static void main(String[] args) throws Exception {
28
29			final String usage = "Usage: java MatrixMultiply <threads> <matrix size (must be a power of two)> [<granularity>] \n Size and granularity must be powers of two.\n For example, try java MatrixMultiply 2 512 16";
30
31			int procs = 0;
32			int n = 2048;
33	dl	1.9	int runs = 32;
34	dl	1.1	try {
35			if (args.length > 0)
36			procs = Integer.parseInt(args[0]);
37			if (args.length > 1)
38			n = Integer.parseInt(args[1]);
39	jsr166	1.2	if (args.length > 2)
40	dl	1.1	granularity = Integer.parseInt(args[2]);
41	jsr166	1.2	if (args.length > 3)
42	dl	1.1	runs = Integer.parseInt(args[2]);
43			}
44	jsr166	1.2
45	dl	1.1	catch (Exception e) {
46			System.out.println(usage);
47			return;
48			}
49	jsr166	1.2
50			if ( ((n & (n - 1)) != 0) \|\|
51	dl	1.1	((granularity & (granularity - 1)) != 0) \|\|
52			granularity < 2) {
53			System.out.println(usage);
54			return;
55			}
56	jsr166	1.2
57	dl	1.9	ForkJoinPool pool = (procs == 0) ? ForkJoinPool.commonPool() :
58	dl	1.1	new ForkJoinPool(procs);
59	jsr166	1.2	System.out.println("procs: " + pool.getParallelism() +
60	dl	1.1	" n: " + n + " granularity: " + granularity +
61			" runs: " + runs);
62	jsr166	1.2
63	dl	1.1	float[][] a = new float[n][n];
64			float[][] b = new float[n][n];
65			float[][] c = new float[n][n];
66	jsr166	1.2
67	dl	1.1	for (int i = 0; i < runs; ++i) {
68			init(a, b, n);
69			long start = System.nanoTime();
70	dl	1.9	new Multiplier(a, 0, 0, b, 0, 0, c, 0, 0, n).invoke();
71	dl	1.1	long time = System.nanoTime() - start;
72			double secs = ((double)time) / NPS;
73			Thread.sleep(100);
74	dl	1.9	System.out.printf("Time: %7.3f ", secs);
75			if ((i & 3) == 3) System.out.println();
76	dl	1.1	// check(c, n);
77			}
78			System.out.println(pool.toString());
79	dl	1.9	if (pool != ForkJoinPool.commonPool())
80			pool.shutdown();
81			Thread.sleep(100);
82	dl	1.1	}
83
84			// To simplify checking, fill with all 1's. Answer should be all n's.
85			static void init(float[][] a, float[][] b, int n) {
86			for (int i = 0; i < n; ++i) {
87			for (int j = 0; j < n; ++j) {
88			a[i][j] = 1.0F;
89			b[i][j] = 1.0F;
90			}
91			}
92			}
93
94			static void check(float[][] c, int n) {
95			for (int i = 0; i < n; i++ ) {
96			for (int j = 0; j < n; j++ ) {
97			if (c[i][j] != n) {
98			throw new Error("Check Failed at [" + i +"]["+j+"]: " + c[i][j]);
99			}
100			}
101			}
102			}
103
104	jsr166	1.2	/**
105	dl	1.1	* Multiply matrices AxB by dividing into quadrants, using algorithm:
106			* <pre>
107	jsr166	1.2	* A x B
108	dl	1.1	*
109	jsr166	1.2	* A11 \| A12 B11 \| B12 A11B11 \| A11B12 A12B21 \| A12B22
110	dl	1.1	* \|----+----\| x \|----+----\| = \|--------+--------\| + \|---------+-------\|
111	jsr166	1.2	* A21 \| A22 B21 \| B21 A21B11 \| A21B21 A22B21 \| A22B22
112	dl	1.1	* </pre>
113			*/
114			static class Multiplier extends RecursiveAction {
115			final float[][] A; // Matrix A
116			final int aRow; // first row of current quadrant of A
117			final int aCol; // first column of current quadrant of A
118
119			final float[][] B; // Similarly for B
120			final int bRow;
121			final int bCol;
122
123			final float[][] C; // Similarly for result matrix C
124			final int cRow;
125			final int cCol;
126
127			final int size; // number of elements in current quadrant
128	jsr166	1.2
129	dl	1.1	Multiplier(float[][] A, int aRow, int aCol,
130			float[][] B, int bRow, int bCol,
131			float[][] C, int cRow, int cCol,
132			int size) {
133			this.A = A; this.aRow = aRow; this.aCol = aCol;
134			this.B = B; this.bRow = bRow; this.bCol = bCol;
135			this.C = C; this.cRow = cRow; this.cCol = cCol;
136			this.size = size;
137			}
138
139			public void compute() {
140
141			if (size <= granularity) {
142			multiplyStride2();
143			}
144
145			else {
146			int h = size / 2;
147
148			invokeAll(new Seq2[] {
149			seq(new Multiplier(A, aRow, aCol, // A11
150			B, bRow, bCol, // B11
151			C, cRow, cCol, // C11
152			h),
153			new Multiplier(A, aRow, aCol+h, // A12
154			B, bRow+h, bCol, // B21
155			C, cRow, cCol, // C11
156			h)),
157	jsr166	1.2
158	dl	1.1	seq(new Multiplier(A, aRow, aCol, // A11
159			B, bRow, bCol+h, // B12
160			C, cRow, cCol+h, // C12
161			h),
162			new Multiplier(A, aRow, aCol+h, // A12
163			B, bRow+h, bCol+h, // B22
164			C, cRow, cCol+h, // C12
165			h)),
166	jsr166	1.2
167	dl	1.1	seq(new Multiplier(A, aRow+h, aCol, // A21
168			B, bRow, bCol, // B11
169			C, cRow+h, cCol, // C21
170			h),
171			new Multiplier(A, aRow+h, aCol+h, // A22
172			B, bRow+h, bCol, // B21
173			C, cRow+h, cCol, // C21
174			h)),
175	jsr166	1.2
176	dl	1.1	seq(new Multiplier(A, aRow+h, aCol, // A21
177			B, bRow, bCol+h, // B12
178			C, cRow+h, cCol+h, // C22
179			h),
180			new Multiplier(A, aRow+h, aCol+h, // A22
181			B, bRow+h, bCol+h, // B22
182			C, cRow+h, cCol+h, // C22
183			h))
184			});
185			}
186			}
187
188	jsr166	1.2	/**
189	dl	1.1	* Version of matrix multiplication that steps 2 rows and columns
190			* at a time. Adapted from Cilk demos.
191			* Note that the results are added into C, not just set into C.
192			* This works well here because Java array elements
193			* are created with all zero values.
194	jsr166	1.3	*/
195	dl	1.1	void multiplyStride2() {
196			for (int j = 0; j < size; j+=2) {
197			for (int i = 0; i < size; i +=2) {
198
199			float[] a0 = A[aRow+i];
200			float[] a1 = A[aRow+i+1];
201	jsr166	1.2
202			float s00 = 0.0F;
203			float s01 = 0.0F;
204			float s10 = 0.0F;
205			float s11 = 0.0F;
206	dl	1.1
207			for (int k = 0; k < size; k+=2) {
208
209			float[] b0 = B[bRow+k];
210
211			s00 += a0[aCol+k] * b0[bCol+j];
212			s10 += a1[aCol+k] * b0[bCol+j];
213			s01 += a0[aCol+k] * b0[bCol+j+1];
214			s11 += a1[aCol+k] * b0[bCol+j+1];
215
216			float[] b1 = B[bRow+k+1];
217
218			s00 += a0[aCol+k+1] * b1[bCol+j];
219			s10 += a1[aCol+k+1] * b1[bCol+j];
220			s01 += a0[aCol+k+1] * b1[bCol+j+1];
221			s11 += a1[aCol+k+1] * b1[bCol+j+1];
222			}
223
224			C[cRow+i] [cCol+j] += s00;
225			C[cRow+i] [cCol+j+1] += s01;
226			C[cRow+i+1][cCol+j] += s10;
227			C[cRow+i+1][cCol+j+1] += s11;
228			}
229			}
230			}
231
232			}
233
234	jsr166	1.2	static Seq2 seq(RecursiveAction task1,
235			RecursiveAction task2) {
236			return new Seq2(task1, task2);
237	dl	1.1	}
238
239			static final class Seq2 extends RecursiveAction {
240			final RecursiveAction fst;
241			final RecursiveAction snd;
242			public Seq2(RecursiveAction task1, RecursiveAction task2) {
243			fst = task1;
244			snd = task2;
245			}
246			public void compute() {
247			fst.invoke();
248			snd.invoke();
249			}
250			}
251			}