Jacobi


1   // Jacobi iteration on a mesh. Based loosely on a Filaments demo
2   
3   import EDU.oswego.cs.dl.util.concurrent.*;
4   
5   public class Jacobi {
6   
7     static final int DEFAULT_LEAFCELLS = 1024;
8   
9     /** 
10     * The maximum number of matrix cells 
11     * at which to stop recursing down and instead directly update.
12     **/
13  
14    static final double EPSILON = 0.001;  // convergence criterion
15  
16    public static void main(String  [] args) {
17      try {
18        int procs;
19        int n;
20        int steps;
21        int granularity = DEFAULT_LEAFCELLS;
22  
23        try {
24          procs = Integer.parseInt(args[0]);
25          n = Integer.parseInt(args[1]);
26          steps = Integer.parseInt(args[2]);
27          if (args.length > 3) granularity = Integer.parseInt(args[3]);
28        }
29  
30        catch (Exception   e) {
31          System.out.println("Usage: java Jacobi <threads> <matrix size> <max steps> [<leafcells>]");
32          return;
33        }
34  
35        // allocate enough space for edges
36  
37        double[][] a = new double[n+2][n+2];
38        double[][] b = new double[n+2][n+2];
39  
40  
41        // Simple initialization for demo. Fill all edges with 1's.
42        // (All interiors are already default-initialized to zero.)
43  
44        for (int k = 0; k < n+2; ++k) {
45          a[k][0] = 1.0;
46          a[k][n+1] = 1.0;
47          a[0][k] = 1.0;
48          a[n+1][k] = 1.0;
49  
50          b[k][0] = 1.0;
51          b[k][n+1] = 1.0;
52          b[0][k] = 1.0;
53          b[n+1][k] = 1.0;
54        }
55  
56        Driver driver = new Driver(a, b, 1, n, 1, n, steps, granularity);
57  
58        FJTaskRunnerGroup g = new FJTaskRunnerGroup(procs);
59        g.invoke(driver);
60        g.stats();
61  
62      }
63      catch (InterruptedException   ex) {}
64    }
65  
66    abstract static class MatrixTree extends FJTask {
67  
68      // maximum difference between old and new values
69      volatile double maxDiff; 
70  
71    }
72  
73  
74    static class LeafNode extends MatrixTree {
75      final double[][] A; // matrix to get old values from
76      final double[][] B; // matrix to put new values into
77  
78      // indices of current submatrix
79      final int loRow;    final int hiRow;
80      final int loCol;    final int hiCol;
81  
82      int steps = 0; // track even/odd steps
83  
84      LeafNode(double[][] A, double[][] B, 
85               int loRow, int hiRow,
86               int loCol, int hiCol) {
87        this.A = A;   this.B = B;
88        this.loRow = loRow; this.hiRow = hiRow;
89        this.loCol = loCol; this.hiCol = hiCol;
90      }
91  
92      public synchronized void run() { 
93        
94        boolean AtoB = (steps++ % 2) == 0;
95        double[][] a = (AtoB)? A : B;
96        double[][] b = (AtoB)? B : A;
97  
98        double md = 0.0; // local for computing max diff
99  
100       for (int i = loRow; i <= hiRow; ++i) {
101         for (int j = loCol; j <= hiCol; ++j) {
102           double v = 0.25 * (a[i-1][j] + a[i][j-1] +
103                              a[i+1][j] + a[i][j+1]);
104           b[i][j] = v;
105 
106           double diff = v - a[i][j];
107           if (diff < 0) diff = -diff;
108           if (diff > md) md = diff;
109         }
110       }
111 
112       maxDiff = md;
113     }
114   }
115 
116 
117 
118   static class FourNode extends MatrixTree {
119     final MatrixTree[] quads;
120 
121     FourNode(MatrixTree q1, MatrixTree q2, 
122              MatrixTree q3, MatrixTree q4) {
123       quads = new MatrixTree[] { q1, q2, q3, q4 };
124     }
125 
126     public void run() { 
127       coInvoke(quads);
128 
129       double md = quads[0].maxDiff;
130       quads[0].reset(); 
131 
132       double m = quads[1].maxDiff;
133       quads[1].reset(); 
134       if (m > md) md = m;
135 
136       m = quads[2].maxDiff;
137       quads[2].reset(); 
138       if (m > md) md = m;
139 
140       m = quads[3].maxDiff;
141       quads[3].reset(); 
142       maxDiff =  (m > md)? m : md;
143     }
144   }
145         
146 
147   static class TwoNode extends MatrixTree {
148     final MatrixTree q1;
149     final MatrixTree q2;
150 
151     TwoNode(MatrixTree q1, MatrixTree q2) {
152       this.q1 = q1; this.q2 = q2;
153     }
154 
155     public void run() { 
156       FJTask.coInvoke(q1, q2);
157       double m1 = q1.maxDiff;
158       double m2 = q2.maxDiff;
159       maxDiff = (m1 > m2)? m1: m2;
160       q1.reset();
161       q2.reset();
162     }
163        
164   }
165         
166 
167   static class Driver extends FJTask {
168     final MatrixTree mat;
169     final int steps;
170 
171     Driver(double[][] A, double[][] B, 
172            int firstRow, int lastRow,
173            int firstCol, int lastCol,
174            int steps, int leafCells) {
175       this.steps = steps;
176       mat = build(A, B, firstRow, lastRow, firstCol, lastCol, leafCells);
177     }
178 
179     MatrixTree build(double[][] a, double[][] b,
180                      int lr, int hr, int lc, int hc, int gran) {
181       int rows = (hr - lr + 1);
182       int cols = (hc - lc + 1);
183 
184       int mr = (lr + hr) / 2; // midpoints
185       int mc = (lc + hc) / 2;
186       
187       int hrows = (mr - lr + 1);
188       int hcols = (mc - lc + 1);
189 
190       if (rows * cols <= gran) {
191         return new LeafNode(a, b, lr, hr, lc, hc);
192       }
193       else if (hrows * hcols >= gran) {
194         return new FourNode(build(a, b, lr,   mr, lc,   mc, gran),
195                             build(a, b, lr,   mr, mc+1, hc, gran),
196                             build(a, b, mr+1, hr, lc,   mc, gran),
197                             build(a, b, mr+1, hr, mc+1, hc, gran));
198       }
199       else if (cols >= rows) {
200         return new TwoNode(build(a, b, lr, hr, lc,   mc, gran),
201                            build(a, b, lr, hr, mc+1, hc, gran));
202       }
203       else {
204         return new TwoNode(build(a, b, lr,   mr, lc, hc, gran),
205                            build(a, b, mr+1, hr, lc, hc, gran));
206         
207       }
208     }
209 
210 
211     public void run() {
212       double md = 0.0;
213 
214       for (int i = 1; i <= steps; ++i) {
215         invoke(mat);
216         md = mat.maxDiff;
217         if (md < EPSILON) { 
218           System.out.println("Converged after " + i + " steps");
219           return;
220         }
221         else
222           mat.reset();
223       }
224 
225       System.out.println("max diff after " + steps + " steps = " + md);
226     }
227   }
228 
229 
230 }
231
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags