KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > org > jscience > mathematics > matrices > 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 JavaDoc;
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 JavaDoc<Operable> NUMERIC_COMPARATOR = new Comparator JavaDoc<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 JavaDoc<Operable>> PIVOT_COMPARATOR
68         = new LocalReference<Comparator JavaDoc<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 JavaDoc("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 JavaDoc<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 JavaDoc<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 }
Popular Tags