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 is7  * 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  *

This class represents the decomposition of a {@link Matrix matrix} into21  * a product of lower and upper triangular matrices22  * (L and U, respectively).

23  * 24  *

This decomposition is typically used to resolve linear systems25  * of equations (Gaussian elimination) or to calculate the determinant26  * of a square {@link Matrix} (O(m³)).

27  * 28  *

Numerical stability is guaranteed through pivoting if the29  * {@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 setting32  * the {@link LocalContext context-local} pivot comparator (see 33  * {@link #setPivotComparator}).

34  * 35  *

Instances of this class are created using the {@link Matrix#lu()}36  * factory method.

37  *38  * @author Jean-Marie Dautelle39  * @version 2.0, June 6, 200440  * @see 41  * LU Decomposition -- from MathWorld 42  */43 public final class LUDecomposition> extends44         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 NUMERIC_COMPARATOR = new Comparator () {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))) // Zero57 return -1;58             if (right.equals(right.plus(right))) // Zero59 return 1;60             return 0;61         }62     };63 64     /**65      * Holds the local comparator.66      */67     private static final LocalReference> PIVOT_COMPARATOR68         = new LocalReference>(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 > LUDecomposition valueOf(Matrix 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 : " + dimension150                     + " 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 cmp) {167         PIVOT_COMPARATOR.set(cmp);168     }169 170     /**171      * Returns the solution X of the equation: A * X = B with172      * this = A.lu() using back and forward substitutions.173      *174      * @param B the input vector.175      * @return the solution X = (1 / A) * B.176      * @throws MatrixException if177      * A.getRowDimension() != B.getRowDimension().178      */179     public Matrix solve(Matrix B) {180         if (_n != B.m)181             throw new MatrixException("Input vector has " + B.m182                     + " rows instead of " + _n);183 184         // Copies B with pivoting.185 final int n = B.n;186         Matrix 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 with224      * this = A.lu() using back and forward substitutions.225      * This method uses {@link PoolContext} internally.226      *227      * @return this.solve(Identity)228      */229     public Matrix inverse() {230         // Calculates inv(U).231 final int n = _n;232         Matrix 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 this291      * {@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 the306      * 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 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 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 largest323 // elements of the matrix source in each column lies324 // 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 important347 // 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 FACTORY_2 = new Factory() {372         protected LUDecomposition create() {373             return new LUDecomposition(1 << 2);374         }375     };376 377     private static final Factory FACTORY_3 = new Factory() {378         protected LUDecomposition create() {379             return new LUDecomposition(1 << 3);380         }381     };382 383     private static final Factory FACTORY_4 = new Factory() {384         protected LUDecomposition create() {385             return new LUDecomposition(1 << 4);386         }387     };388 389     private static final Factory FACTORY_5 = new Factory() {390         protected LUDecomposition create() {391             return new LUDecomposition(1 << 5);392         }393     };394 395     private static final Factory FACTORY_6 = new Factory() {396         protected LUDecomposition create() {397             return new LUDecomposition(1 << 6);398         }399     };400 401     private static final Factory FACTORY_7 = new Factory() {402         protected LUDecomposition create() {403             return new LUDecomposition(1 << 7);404         }405     };406 407     private static final Factory FACTORY_8 = new Factory() {408         protected LUDecomposition create() {409             return new LUDecomposition(1 << 8);410         }411     };412 413     private static final Factory FACTORY_9 = new Factory() {414         protected LUDecomposition create() {415             return new LUDecomposition(1 << 9);416         }417     };418 419     private static final Factory FACTORY_10 = new Factory() {420         protected LUDecomposition create() {421             return new LUDecomposition(1 << 10);422         }423     };424 425     private static final Factory FACTORY_11 = new Factory() {426         protected LUDecomposition create() {427             return new LUDecomposition(1 << 11);428         }429     };430 431     private static final Factory FACTORY_12 = new Factory() {432         protected LUDecomposition create() {433             return new LUDecomposition(1 << 12);434         }435     };436 437     private static final Factory FACTORY_13 = new Factory() {438         protected LUDecomposition create() {439             return new LUDecomposition(1 << 13);440         }441     };442 443     private static final Factory FACTORY_14 = new Factory() {444         protected LUDecomposition create() {445             return new LUDecomposition(1 << 14);446         }447     };448 449     private static final Factory FACTORY_15 = new Factory() {450         protected LUDecomposition create() {451             return new LUDecomposition(1 << 15);452         }453     };454 455 }` Popular Tags