|
1
2 import EDU.oswego.cs.dl.util.concurrent.*;
3
4
8
9 public class LU {
10
11 static final int BLOCK_SIZE = 16;
13
14 static final boolean CHECK = false;
16 public static void main(String  [] args) {
17
18 final String usage = "Usage: java LU <threads> <matrix size (must be a power of two)> [runs] \n For example, try java LU 2 512";
19
20 try {
21 int procs;
22 int n;
23 int runs = 1;
24 try {
25 procs = Integer.parseInt(args[0]);
26 n = Integer.parseInt(args[1]);
27 if (args.length > 2) runs = Integer.parseInt(args[2]);
28
29 }
30
31 catch (Exception e) {
32 System.out.println(usage);
33 return;
34 }
35
36 if ( ((n & (n - 1)) != 0)) {
37 System.out.println(usage);
38 return;
39 }
40
41 for (int run = 0; run < runs; ++run) {
42
43 double[][] m = new double[n][n];
44 randomInit(m, n);
45
46
47 double[][] copy = null;
48 if (CHECK) {
49 copy = new double[n][n];
50 for (int i = 0; i < n; ++i) {
51 for (int j = 0; j < n; ++j) {
52 copy[i][j] = m[i][j];
53 }
54 }
55 }
56
57 Block M = new Block(m, 0, 0);
58
59 FJTaskRunnerGroup g = new FJTaskRunnerGroup(procs);
60 g.invoke(new LowerUpper(n, M));
61 g.stats();
62 g.interruptAll();
63 if (CHECK) check(m, copy, n);
64 }
65 }
66 catch (InterruptedException ex) {}
67 }
68
69
70 static void randomInit(double[][] M, int n) {
71
72 java.util.Random rng = new java.util.Random ();
73
74 for (int i = 0; i < n; ++i)
75 for (int j = 0; j < n; ++j)
76 M[i][j] = rng.nextDouble();
77
78 for (int k = 0; k < n; ++k)
80 M[k][k] *= 10.0;
81 }
82
83 static void check(double[][] LU, double[][] M, int n) {
84
85 double maxDiff = 0.0;
87 for (int i = 0; i < n; ++i) {
88 for (int j = 0; j < n; ++j) {
89 double v = 0.0;
90 int k;
91 for (k = 0; k < i && k <= j; k++ ) v += LU[i][k] * LU[k][j];
92 if (k == i && k <= j ) v += LU[k][j];
93 double diff = M[i][j] - v;
94 if (diff < 0) diff = -diff;
95 if (diff > 0.001) {
96 System.out.println("large diff at[" + i + "," + j + "]: " + M[i][j] + " vs " + v);
97 }
98 if (diff > maxDiff) maxDiff = diff;
99 }
100 }
101
102 System.out.println("Max difference = " + maxDiff);
103 }
104
105
106 static class Block {
108 final double[][] m;
109 final int loRow;
110 final int loCol;
111
112 Block(double[][] mat, int lr, int lc) {
113 m = mat; loRow = lr; loCol = lc;
114 }
115 }
116
117 static class Schur extends FJTask {
118 final int size;
119 final Block V;
120 final Block W;
121 final Block M;
122
123 Schur(int size, Block V, Block W, Block M) {
124 this.size = size; this.V = V; this.W = W; this.M = M;
125 }
126
127 void schur() { for (int j = 0; j < BLOCK_SIZE; ++j) {
129 for (int i = 0; i < BLOCK_SIZE; ++i) {
130 double s = M.m[i+M.loRow][j+M.loCol];
131 for (int k = 0; k < BLOCK_SIZE; ++k) {
132 s -= V.m[i+V.loRow][k+V.loCol] * W.m[k+W.loRow][j+W.loCol];
133 }
134 M.m[i+M.loRow][j+M.loCol] = s;
135 }
136 }
137 }
138
139 public void run() {
140 if (size == BLOCK_SIZE) {
141 schur();
142 }
143 else {
144 int h = size / 2;
145
146 Block M00 = new Block(M.m, M.loRow, M.loCol);
147 Block M01 = new Block(M.m, M.loRow, M.loCol+h);
148 Block M10 = new Block(M.m, M.loRow+h, M.loCol);
149 Block M11 = new Block(M.m, M.loRow+h, M.loCol+h);
150
151 Block V00 = new Block(V.m, V.loRow, V.loCol);
152 Block V01 = new Block(V.m, V.loRow, V.loCol+h);
153 Block V10 = new Block(V.m, V.loRow+h, V.loCol);
154 Block V11 = new Block(V.m, V.loRow+h, V.loCol+h);
155
156 Block W00 = new Block(W.m, W.loRow, W.loCol);
157 Block W01 = new Block(W.m, W.loRow, W.loCol+h);
158 Block W10 = new Block(W.m, W.loRow+h, W.loCol);
159 Block W11 = new Block(W.m, W.loRow+h, W.loCol+h);
160
161 coInvoke(new FJTask[] {
162 seq(new Schur(h, V00, W00, M00),
163 new Schur(h, V01, W10, M00)),
164 seq(new Schur(h, V00, W01, M01),
165 new Schur(h, V01, W11, M01)),
166 seq(new Schur(h, V10, W00, M10),
167 new Schur(h, V11, W10, M10)),
168 seq(new Schur(h, V10, W01, M11),
169 new Schur(h, V11, W11, M11))
170 });
171
172 }
173 }
174 }
175
176 static class Lower extends FJTask {
177
178 final int size;
179 final Block L;
180 final Block M;
181
182 Lower(int size, Block L, Block M) {
183 this.size = size; this.L = L; this.M = M;
184 }
185
186
187 void lower() { for (int i = 1; i < BLOCK_SIZE; ++i) {
189 for (int k = 0; k < i; ++k) {
190 double a = L.m[i+L.loRow][k+L.loCol];
191 double[] x = M.m[k+M.loRow];
192 double[] y = M.m[i+M.loRow];
193 int n = BLOCK_SIZE;
194 for (int p = n-1; p >= 0; --p) {
195 y[p+M.loCol] -= a * x[p+M.loCol];
196 }
197 }
198 }
199 }
200
201 public void run() {
202 if (size == BLOCK_SIZE) {
203 lower();
204 }
205 else {
206 int h = size / 2;
207
208 Block M00 = new Block(M.m, M.loRow, M.loCol);
209 Block M01 = new Block(M.m, M.loRow, M.loCol+h);
210 Block M10 = new Block(M.m, M.loRow+h, M.loCol);
211 Block M11 = new Block(M.m, M.loRow+h, M.loCol+h);
212
213 Block L00 = new Block(L.m, L.loRow, L.loCol);
214 Block L01 = new Block(L.m, L.loRow, L.loCol+h);
215 Block L10 = new Block(L.m, L.loRow+h, L.loCol);
216 Block L11 = new Block(L.m, L.loRow+h, L.loCol+h);
217
218 coInvoke(
219 new Seq(new FJTask[] {
220 new Lower(h, L00, M00),
221 new Schur(h, L10, M00, M10),
222 new Lower(h, L11, M10)
223 }),
224 new Seq(new FJTask[] {
225 new Lower(h, L00, M01),
226 new Schur(h, L10, M01, M11),
227 new Lower(h, L11, M11)
228 })
229 );
230 }
231 }
232 }
233
234
235 static class Upper extends FJTask {
236
237 final int size;
238 final Block U;
239 final Block M;
240
241 Upper(int size, Block U, Block M) {
242 this.size = size; this.U = U; this.M = M;
243 }
244
245 void upper() { for (int i = 0; i < BLOCK_SIZE; ++i) {
247 for (int k = 0; k < BLOCK_SIZE; ++k) {
248 double a = M.m[i+M.loRow][k+M.loCol] / U.m[k+U.loRow][k+U.loCol];
249 M.m[i+M.loRow][k+M.loCol] = a;
250 double[] x = U.m[k+U.loRow];
251 double[] y = M.m[i+M.loRow];
252 int n = BLOCK_SIZE - k - 1;
253 for (int p = n - 1; p >= 0; --p) {
254 y[p+k+1+M.loCol] -= a * x[p+k+1+U.loCol];
255 }
256 }
257 }
258 }
259
260
261 public void run() {
262 if (size == BLOCK_SIZE) {
263 upper();
264 }
265 else {
266 int h = size / 2;
267
268 Block M00 = new Block(M.m, M.loRow, M.loCol);
269 Block M01 = new Block(M.m, M.loRow, M.loCol+h);
270 Block M10 = new Block(M.m, M.loRow+h, M.loCol);
271 Block M11 = new Block(M.m, M.loRow+h, M.loCol+h);
272
273 Block U00 = new Block(U.m, U.loRow, U.loCol);
274 Block U01 = new Block(U.m, U.loRow, U.loCol+h);
275 Block U10 = new Block(U.m, U.loRow+h, U.loCol);
276 Block U11 = new Block(U.m, U.loRow+h, U.loCol+h);
277
278 coInvoke(
279 new Seq(new FJTask[] {
280 new Upper(h, U00, M00),
281 new Schur(h, M00, U01, M01),
282 new Upper(h, U11, M01)
283 }),
284 new Seq(new FJTask[] {
285 new Upper(h, U00, M10),
286 new Schur(h, M10, U01, M11),
287 new Upper(h, U11, M11)
288 })
289 );
290 }
291 }
292 }
293
294
295 static class LowerUpper extends FJTask {
296
297 final int size;
298 final Block M;
299
300 LowerUpper(int size, Block M) {
301 this.size = size; this.M = M;
302 }
303
304 void lu() { for (int k = 0; k < BLOCK_SIZE; ++k) {
306 for (int i = k+1; i < BLOCK_SIZE; ++i) {
307 double b = M.m[k+M.loRow][k+M.loCol];
308 double a = M.m[i+M.loRow][k+M.loCol] / b;
309 M.m[i+M.loRow][k+M.loCol] = a;
310 double[] x = M.m[k+M.loRow];
311 double[] y = M.m[i+M.loRow];
312 int n = BLOCK_SIZE-k-1;
313 for (int p = n-1; p >= 0; --p) {
314 y[k+1+p+M.loCol] -= a * x[k+1+p+M.loCol];
315 }
316 }
317 }
318 }
319
320 public void run() {
321 if (size == BLOCK_SIZE) {
322 lu();
323 }
324 else {
325 int h = size / 2;
326
327 Block M00 = new Block(M.m, M.loRow, M.loCol);
328 Block M01 = new Block(M.m, M.loRow, M.loCol+h);
329 Block M10 = new Block(M.m, M.loRow+h, M.loCol);
330 Block M11 = new Block(M.m, M.loRow+h, M.loCol+h);
331
332
333 invoke(new LowerUpper(h, M00));
334
335 coInvoke(new Lower(h, M00, M01),
336 new Upper(h, M00, M10));
337
338 invoke(new Schur(h, M10, M01, M11));
339
340 invoke(new LowerUpper(h, M11));
341
342 }
343 }
344 }
345 }
346
347 |
Popular Tags |
Java API By Example, From Geeks To Geeks.