LUDecomposition


1   /*
2    * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
3    * Copyright (C) 2005 - JScience (http://jscience.org/)
4    * All rights reserved.
5    * 
6    * Permission to use, copy, modify, and distribute this software is
7    * freely granted, provided that this notice is preserved.
8    */
9   package org.jscience.mathematics.matrices;
10  
11  import org.jscience.mathematics.numbers.Numeric;
12  import java.util.Comparator  ;
13  
14  import javolution.realtime.LocalContext;
15  import javolution.realtime.LocalReference;
16  import javolution.realtime.PoolContext;
17  import javolution.realtime.RealtimeObject;
18  
19  /**
20   * <p> This class represents the decomposition of a {@link Matrix matrix} into
21   *     a product of lower and upper triangular matrices
22   *     (L and U, respectively).</p>
23   *     
24   * <p> This decomposition</a> is typically used to resolve linear systems
25   *     of equations (Gaussian elimination) or to calculate the determinant
26   *     of a square {@link Matrix} (<code>O(m³)</code>).</p>
27   *     
28   * <p> Numerical stability is guaranteed through pivoting if the
29   *     {@link Operable} elements of this matrix are 
30   *     {@link org.jscience.mathematics.numbers.Numeric Numeric}
31   *     For others elements types, numerical stability can be ensured by setting
32   *     the {@link LocalContext context-local} pivot comparator (see 
33   *     {@link #setPivotComparator}).</p>
34   *     
35   * <p> Instances of this class are created using the {@link Matrix#lu()}
36   *     factory method.</p>
37   *
38   * @author  <a HREF="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
39   * @version 2.0, June 6, 2004
40   * @see <a HREF="http://mathworld.wolfram.com/LUDecomposition.html">
41   *      LU Decomposition -- from MathWorld</a> 
42   */
43  public final class LUDecomposition<O extends Operable<O>> extends
44          RealtimeObject {
45  
46      /**
47       * Holds the comparator for pivoting 
48       * {@link org.jscience.mathematics.numbers.Numeric Numeric}
49       * instances.
50       */
51      public static final Comparator  <Operable> NUMERIC_COMPARATOR = new Comparator  <Operable>() {
52  
53          public int compare(Operable left, Operable right) {
54              if ((left instanceof Numeric) && (right instanceof Numeric))
55                  return ((Numeric) left).isLargerThan((Numeric) right) ? 1 : -1;
56              if (left.equals(left.plus(left))) // Zero
57                  return -1;
58              if (right.equals(right.plus(right))) // Zero
59                  return 1;
60              return 0;
61          }
62      };
63  
64      /**
65       * Holds the local comparator.
66       */
67      private static final LocalReference<Comparator  <Operable>> PIVOT_COMPARATOR
68          = new LocalReference<Comparator  <Operable>>(
69              NUMERIC_COMPARATOR);
70  
71      /**
72       * Holds the dimension of the square matrix source.
73       */
74      private int _n;
75  
76      /**
77       * Holds the LU elements.
78       */
79      private O[] _lu;
80  
81      /**
82       * Holds temporary buffer.
83       */
84      private O[] _tmp;
85  
86      /**
87       * Holds the pivots indexes.
88       */
89      private final int[] _pivots;
90  
91      /**
92       * Holds the number of permutation performed.
93       */
94      private int _permutationCount;
95  
96      /**
97       * Base constructor.
98       * 
99       * @param dimension the dimension of the square matrix source.
100      */
101     private LUDecomposition(int dimension) {
102         _lu = (O[]) new Operable[dimension * dimension];
103         _tmp = (O[]) new Operable[dimension * dimension];
104         _pivots = new int[dimension];
105     }
106 
107     /**
108      * Returns the LU decomposition of the specified matrix..
109      *
110      * @param  source the matrix for which the LU decomposition is 
111      *         returned.
112      * @return the corresponding LU decomposition.
113      * @throws MatrixException if this matrix is not square.
114      */
115     static <O extends Operable<O>> LUDecomposition<O> valueOf(Matrix<O> source) {
116         if (!source.isSquare())
117             throw new MatrixException("Matrix is not square");
118         int dimension = source.m;
119         LUDecomposition lu;
120         if (dimension <= 1 << 2) {
121             lu = FACTORY_2.object();
122         } else if (dimension <= 1 << 3) {
123             lu = FACTORY_3.object();
124         } else if (dimension <= 1 << 4) {
125             lu = FACTORY_4.object();
126         } else if (dimension <= 1 << 5) {
127             lu = FACTORY_5.object();
128         } else if (dimension <= 1 << 6) {
129             lu = FACTORY_6.object();
130         } else if (dimension <= 1 << 7) {
131             lu = FACTORY_7.object();
132         } else if (dimension <= 1 << 8) {
133             lu = FACTORY_8.object();
134         } else if (dimension <= 1 << 9) {
135             lu = FACTORY_9.object();
136         } else if (dimension <= 1 << 10) {
137             lu = FACTORY_10.object();
138         } else if (dimension <= 1 << 11) {
139             lu = FACTORY_11.object();
140         } else if (dimension <= 1 << 12) {
141             lu = FACTORY_12.object();
142         } else if (dimension <= 1 << 13) {
143             lu = FACTORY_13.object();
144         } else if (dimension <= 1 << 14) {
145             lu = FACTORY_14.object();
146         } else if (dimension <= 1 << 15) {
147             lu = FACTORY_15.object();
148         } else {
149             throw new UnsupportedOperationException  ("Dimension : " + dimension
150                     + " too large");
151         }
152         lu._n = dimension;
153         lu._permutationCount = 0;
154         lu.construct(source);
155         return lu;
156 
157     }
158 
159     /**
160      * Sets the {@link LocalContext local} comparator used 
161      * for pivoting during LU decomposition (default 
162      * {@link #NUMERIC_COMPARATOR}).
163      *
164      * @param  cmp the comparator for pivoting.
165      */
166     public void setPivotComparator(Comparator  <Operable> cmp) {
167         PIVOT_COMPARATOR.set(cmp);
168     }
169 
170     /**
171      * Returns the solution X of the equation: A * X = B  with
172      * <code>this = A.lu()</code> using back and forward substitutions.
173      *
174      * @param  B the input vector.
175      * @return the solution X = (1 / A) * B.
176      * @throws MatrixException if
177      *         <code>A.getRowDimension() != B.getRowDimension()</code>.
178      */
179     public Matrix<O> solve(Matrix<O> B) {
180         if (_n != B.m)
181             throw new MatrixException("Input vector has " + B.m
182                     + " rows instead of " + _n);
183 
184         // Copies B with pivoting.
185         final int n = B.n;
186         Matrix<O> X = Matrix.newInstance(_n, n);
187         final O[] x = X.o;
188         for (int i = 0; i < _n; i++) {
189             for (int j = 0; j < n; j++) {
190                 x[i * n + j] = B.o[_pivots[i] * n + j];
191             }
192         }
193 
194         // Solves L * Y = pivot(B)
195         for (int k = 0; k < _n; k++) {
196             for (int i = k + 1; i < _n; i++) {
197                 O luik = _lu[i * _n + k];
198                 for (int j = 0; j < n; j++) {
199                     x[i * n + j] = x[i * n + j].plus(luik.times(x[k * n + j])
200                             .opposite());
201                 }
202             }
203         }
204 
205         // Solves U * X = Y;
206         for (int k = _n - 1; k >= 0; k--) {
207             for (int j = 0; j < n; j++) {
208                 x[k * n + j] = (_lu[k * _n + k].reciprocal())
209                         .times(x[k * n + j]);
210             }
211             for (int i = 0; i < k; i++) {
212                 O luik = _lu[i * _n + k];
213                 for (int j = 0; j < n; j++) {
214                     x[i * n + j] = x[i * n + j].plus(luik.times(x[k * n + j])
215                             .opposite());
216                 }
217             }
218         }
219         return X;
220     }
221 
222     /**
223      * Returns the solution X of the equation: A * X = Identity  with
224      * <code>this = A.lu()</code> using back and forward substitutions.
225      * This method uses {@link PoolContext} internally.
226      *
227      * @return <code>this.solve(Identity)</code>
228      */
229     public Matrix<O> inverse() {
230         // Calculates inv(U).
231         final int n = _n;
232         Matrix<O> R = Matrix.newInstance(n, n);
233         final O[] r = R.o;
234         for (int i = 0; i < n; i++) {
235             for (int j = i; j < n; j++) {
236                 r[i * n + j] = _lu[i * n + j];
237             }
238         }
239         for (int j = n - 1; j >= 0; j--) {
240             r[j * n + j] = r[j * n + j].reciprocal();
241             for (int i = j - 1; i >= 0; i--) {
242                 //PoolContext.enter();
243                 //try {
244                     O sum = r[i * n + j].times(r[j * n + j]).opposite();
245                     for (int k = j - 1; k > i; k--) {
246                         sum = sum.plus(r[i * n + k].times(r[k * n + j])
247                                 .opposite());
248                     }
249                     r[i * n + j] = (r[i * n + i].reciprocal())
250                             .times(sum);
251                     r[i * n + j].move(ObjectSpace.OUTER);
252                 //} finally {
253                 //    PoolContext.exit();
254                 //}
255             }
256         }
257         // Solves inv(A) * L = inv(U)
258         for (int i = 0; i < n; i++) {
259             for (int j = n - 2; j >= 0; j--) {
260                 //PoolContext.enter();
261                 //try {
262                     for (int k = j + 1; k < n; k++) {
263                         O lukj = _lu[k * n + j];
264                         if (r[i * n + j] != null) {
265                             r[i * n + j] = r[i * n + j].plus(r[i * n + k]
266                                     .times(lukj).opposite());
267                         } else {
268                             r[i * n + j] = r[i * n + k].times(lukj).opposite();
269                         }
270                     }
271                     r[i * n + j].move(ObjectSpace.OUTER);
272                 //} finally {
273                 //    PoolContext.exit();
274                 //}
275             }
276         }
277         // Swaps columns (reverses pivots permutations).
278         for (int i = 0; i < n; i++) {
279             for (int j = 0; j < n; j++) {
280                 _tmp[j] = r[i * n + j];
281             }
282             for (int j = 0; j < n; j++) {
283                 r[i * n + _pivots[j]] = _tmp[j];
284             }
285         }
286         return R;
287     }
288 
289     /**
290      * Returns the determinant of the {@link Matrix} having this
291      * {@link LUDecomposition}.
292      *
293      * @return the determinant of the matrix source.
294      */
295     public O determinant() {
296         O product = _lu[0];
297         for (int i = 1; i < _n; i++) {
298             product = product.times(_lu[i * _n + i]);
299         }
300         return ((_permutationCount & 1) == 0) ? product : product.opposite();
301     }
302 
303     /**
304      * Constructs the LU decomposition of the specified matrix.
305      * We make the choise of Lii = ONE (diagonal elements of the
306      * lower triangular matrix are multiplicative identities).
307      *
308      * @param  source the matrix to decompose.
309      * @throws MatrixException if the matrix source is not square.
310      */
311     private void construct(Matrix<O> source) {
312         System.arraycopy(source.o, 0, _lu, 0, _n * _n);
313         for (int i = 0; i < _n; i++) {
314             _pivots[i] = i;
315         }
316 
317         // Main loop.
318         Comparator  <Operable> cmp = PIVOT_COMPARATOR.get();
319         final int n = _n;
320         for (int k = 0; k < _n; k++) {
321 
322             // Rearranges the rows so that the absolutely largest
323             // elements of the matrix source in each column lies
324             // in the diagonal.
325             int pivot = k;
326             for (int i = k + 1; i < n; i++) {
327                 if (cmp.compare(_lu[i * n + k], _lu[pivot * n + k]) > 0) {
328                     pivot = i;
329                 }
330             }
331             if (pivot != k) { // Exchanges.
332                 for (int j = 0; j < n; j++) {
333                     O tmp = _lu[pivot * n + j];
334                     _lu[pivot * n + j] = _lu[k * n + j];
335                     _lu[k * n + j] = tmp;
336                 }
337                 int j = _pivots[pivot];
338                 _pivots[pivot] = _pivots[k];
339                 _pivots[k] = j;
340                 _permutationCount++;
341             }
342 
343             // Computes multipliers and eliminate k-th column.
344             O lukkInv = _lu[k * n + k].reciprocal();
345             for (int i = k + 1; i < n; i++) {
346                 // Multiplicative order is important
347                 // for non-commutative elements.
348                 _lu[i * n + k] = _lu[i * n + k].times(lukkInv);
349                 for (int j = k + 1; j < n; j++) {
350                     _lu[i * n + j] = _lu[i * n + j].plus(_lu[i * n + k].times(
351                             _lu[k * n + j]).opposite());
352                 }
353             }
354         }
355     }
356 
357     // Overrides.
358     public boolean move(ObjectSpace os) {
359         if (super.move(os)) {
360             for (int i = _n * _n; i > 0;) {
361                 _lu[--i].move(os);
362             }
363             return true;
364         }
365         return false;
366     }
367 
368     // TBD: Use recursive structures instead of large arrays.
369     //
370 
371     private static final Factory<LUDecomposition> FACTORY_2 = new Factory<LUDecomposition>() {
372         protected LUDecomposition create() {
373             return new LUDecomposition(1 << 2);
374         }
375     };
376 
377     private static final Factory<LUDecomposition> FACTORY_3 = new Factory<LUDecomposition>() {
378         protected LUDecomposition create() {
379             return new LUDecomposition(1 << 3);
380         }
381     };
382 
383     private static final Factory<LUDecomposition> FACTORY_4 = new Factory<LUDecomposition>() {
384         protected LUDecomposition create() {
385             return new LUDecomposition(1 << 4);
386         }
387     };
388 
389     private static final Factory<LUDecomposition> FACTORY_5 = new Factory<LUDecomposition>() {
390         protected LUDecomposition create() {
391             return new LUDecomposition(1 << 5);
392         }
393     };
394 
395     private static final Factory<LUDecomposition> FACTORY_6 = new Factory<LUDecomposition>() {
396         protected LUDecomposition create() {
397             return new LUDecomposition(1 << 6);
398         }
399     };
400 
401     private static final Factory<LUDecomposition> FACTORY_7 = new Factory<LUDecomposition>() {
402         protected LUDecomposition create() {
403             return new LUDecomposition(1 << 7);
404         }
405     };
406 
407     private static final Factory<LUDecomposition> FACTORY_8 = new Factory<LUDecomposition>() {
408         protected LUDecomposition create() {
409             return new LUDecomposition(1 << 8);
410         }
411     };
412 
413     private static final Factory<LUDecomposition> FACTORY_9 = new Factory<LUDecomposition>() {
414         protected LUDecomposition create() {
415             return new LUDecomposition(1 << 9);
416         }
417     };
418 
419     private static final Factory<LUDecomposition> FACTORY_10 = new Factory<LUDecomposition>() {
420         protected LUDecomposition create() {
421             return new LUDecomposition(1 << 10);
422         }
423     };
424 
425     private static final Factory<LUDecomposition> FACTORY_11 = new Factory<LUDecomposition>() {
426         protected LUDecomposition create() {
427             return new LUDecomposition(1 << 11);
428         }
429     };
430 
431     private static final Factory<LUDecomposition> FACTORY_12 = new Factory<LUDecomposition>() {
432         protected LUDecomposition create() {
433             return new LUDecomposition(1 << 12);
434         }
435     };
436 
437     private static final Factory<LUDecomposition> FACTORY_13 = new Factory<LUDecomposition>() {
438         protected LUDecomposition create() {
439             return new LUDecomposition(1 << 13);
440         }
441     };
442 
443     private static final Factory<LUDecomposition> FACTORY_14 = new Factory<LUDecomposition>() {
444         protected LUDecomposition create() {
445             return new LUDecomposition(1 << 14);
446         }
447     };
448 
449     private static final Factory<LUDecomposition> FACTORY_15 = new Factory<LUDecomposition>() {
450         protected LUDecomposition create() {
451             return new LUDecomposition(1 << 15);
452         }
453     };
454 
455 }
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags