001: /*
002: Copyright © 2007 Stefano Chizzolini. http://clown.stefanochizzolini.it
003:
004: Contributors:
005: * Stefano Chizzolini (original code developer, http://www.stefanochizzolini.it):
006: contributed code is Copyright © 2007 by Stefano Chizzolini.
007:
008: This file should be part of the source code distribution of "PDF Clown library"
009: (the Program): see the accompanying README files for more info.
010:
011: This Program is free software; you can redistribute it and/or modify it under
012: the terms of the GNU General Public License as published by the Free Software
013: Foundation; either version 2 of the License, or (at your option) any later version.
014:
015: This Program is distributed in the hope that it will be useful, but WITHOUT ANY
016: WARRANTY, either expressed or implied; without even the implied warranty of
017: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the License for more details.
018:
019: You should have received a copy of the GNU General Public License along with this
020: Program (see README files); if not, go to the GNU website (http://www.gnu.org/).
021:
022: Redistribution and use, with or without modification, are permitted provided that such
023: redistributions retain the above copyright notice, license and disclaimer, along with
024: this list of conditions.
025: */
026:
027: package it.stefanochizzolini.clown.util.math;
028:
029: /**
030: LU matrix decomposition.
031: <p>The LU decomposition is a lower triangular matrix L, an upper triangular matrix U,
032: and a permutation <code>size</code>-long pivot vector.</p>
033: <h3>Remarks</h3>
034: <p>This class is a specialized adaptation from the original JAMA (Java Matrix) project,
035: brought to the public domain by The MathWorks, Inc. and the National Institute of Standards
036: and Technology (see <a href="http://math.nist.gov/javanumerics/jama/">JAMA home page</a>).</p>
037:
038: @author Stefano Chizzolini
039: @version 0.0.4, 07/12/07
040: @since 0.0.4
041: */
042: public final class LUDecomposition {
043: // <class>
044: // <dynamic>
045: /** Array for internal storage of decomposition. */
046: private double[][] data;
047:
048: /** Matrix size. */
049: private int size;
050:
051: /** Pivot sign. */
052: private int pivsign;
053:
054: /** Internal storage of pivot vector. */
055: private int[] piv;
056:
057: // <constructors>
058: /**
059: @param matrix Matrix to decompose.
060: */
061: public LUDecomposition(SquareMatrix matrix) {
062: /*
063: NOTE: Use a "left-looking", dot-product, Crout/Doolittle algorithm.
064: */
065: this .data = ((SquareMatrix) matrix.clone()).getData();
066: this .size = matrix.getSize();
067:
068: this .piv = new int[size];
069: for (int i = 0; i < size; i++) {
070: piv[i] = i;
071: }
072:
073: pivsign = 1;
074: double[] LUrowi;
075: double[] LUcolj = new double[size];
076:
077: for (int j = 0; j < size; j++) {
078: // Making a copy of the j-th column to localize references...
079: for (int i = 0; i < size; i++) {
080: LUcolj[i] = data[i][j];
081: }
082:
083: // Applying previous transformations...
084: for (int i = 0; i < size; i++) {
085: LUrowi = data[i];
086:
087: int kmax = Math.min(i, j);
088: double s = 0.0;
089: for (int k = 0; k < kmax; k++) {
090: s += LUrowi[k] * LUcolj[k];
091: }
092:
093: LUrowi[j] = LUcolj[i] -= s;
094: }
095:
096: // Finding pivot and exchanging if necessary...
097: int p = j;
098: for (int i = j + 1; i < size; i++) {
099: if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
100: p = i;
101: }
102: }
103:
104: if (p != j) {
105: for (int k = 0; k < size; k++) {
106: double t = data[p][k];
107: data[p][k] = data[j][k];
108: data[j][k] = t;
109: }
110:
111: int k = piv[p];
112: piv[p] = piv[j];
113: piv[j] = k;
114:
115: pivsign = -pivsign;
116: }
117:
118: // Computing multipliers...
119: if (j < size & data[j][j] != 0.0) {
120: for (int i = j + 1; i < size; i++) {
121: data[i][j] /= data[j][j];
122: }
123: }
124: }
125: }
126:
127: // </constructors>
128:
129: /**
130: Gets the determinant.
131: */
132: public double getDet() {
133: double d = (double) pivsign;
134: for (int j = 0; j < size; j++) {
135: d *= data[j][j];
136: }
137:
138: return d;
139: }
140:
141: /**
142: Gets the lower triangular factor.
143: */
144: public SquareMatrix getL() {
145: SquareMatrix X = new SquareMatrix(size);
146: double[][] L = X.getData();
147: for (int i = 0; i < size; i++) {
148: for (int j = 0; j < size; j++) {
149: if (i > j) {
150: L[i][j] = data[i][j];
151: } else if (i == j) {
152: L[i][j] = 1.0;
153: } else {
154: L[i][j] = 0.0;
155: }
156: }
157: }
158:
159: return X;
160: }
161:
162: /**
163: Gets the upper triangular factor.
164: */
165: public SquareMatrix getU() {
166: SquareMatrix X = new SquareMatrix(size);
167: double[][] U = X.getData();
168: for (int i = 0; i < size; i++) {
169: for (int j = 0; j < size; j++) {
170: if (i <= j) {
171: U[i][j] = data[i][j];
172: } else {
173: U[i][j] = 0.0;
174: }
175: }
176: }
177:
178: return X;
179: }
180:
181: /**
182: Gets the pivot permutation vector.
183: */
184: public int[] getPivot() {
185: int[] pivot = new int[size];
186: for (int i = 0; i < size; i++) {
187: pivot[i] = piv[i];
188: }
189:
190: return pivot;
191: }
192:
193: /**
194: Gets whether the matrix is non-singular.
195: */
196: public boolean isNonsingular() {
197: for (int j = 0; j < size; j++)
198: if (data[j][j] == 0)
199: return false;
200:
201: return true;
202: }
203:
204: /**
205: Solve [this] * [return] = [target]
206: @param target Resulting matrix.
207: @return Solution.
208: @exception IllegalArgumentException Matrix sizes must agree.
209: @exception RuntimeException Matrix is singular.
210: */
211: public SquareMatrix solve(SquareMatrix target) {
212: if (target.getSize() != size)
213: throw new IllegalArgumentException(
214: "Matrix size must agree.");
215: if (!this .isNonsingular())
216: throw new RuntimeException("Matrix is singular.");
217:
218: // Copy right hand side with pivoting
219: SquareMatrix Xmat = target.getMatrix(piv, 0);
220: double[][] X = Xmat.getData();
221:
222: // Solve L*Y = target(piv,:)
223: for (int k = 0; k < size; k++) {
224: for (int i = k + 1; i < size; i++) {
225: for (int j = 0; j < size; j++) {
226: X[i][j] -= X[k][j] * data[i][k];
227: }
228: }
229: }
230:
231: // Solve U*X = Y;
232: for (int k = size - 1; k >= 0; k--) {
233: for (int j = 0; j < size; j++) {
234: X[k][j] /= data[k][k];
235: }
236:
237: for (int i = 0; i < k; i++) {
238: for (int j = 0; j < size; j++) {
239: X[i][j] -= X[k][j] * data[i][k];
240: }
241: }
242: }
243:
244: return Xmat;
245: }
246: // </dynamic>
247: // </class>
248: }
|