|
/Users/lyon/j4p/src/math/linearAlgebra/LUDecomposition.java
|
1 package math.linearAlgebra;
2
3 /** LU Decomposition.
4 <P>
5 For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
6 unit lower triangular matrix L, an n-by-n upper triangular matrix U,
7 and a permutation vector piv of length m so that A(piv,:) = L*U.
8 If m < n, then L is m-by-m and U is m-by-n.
9 <P>
10 The LU decompostion with pivoting always exists, even if the matrix is
11 singular, so the constructor will never fail. The primary use of the
12 LU decomposition is in the solution of square systems of simultaneous
13 linear equations. This will fail if isNonsingular() returns false.
14 */
15
16 public class LUDecomposition implements java.io.Serializable {
17
18 /* ------------------------
19 Class variables
20 * ------------------------ */
21
22 /** Array for internal storage of decomposition.
23 @serial internal array storage.
24 */
25 private double[][] LU;
26
27 /** Row and column dimensions, and pivot sign.
28 @serial column dimension.
29 @serial row dimension.
30 @serial pivot sign.
31 */
32 private int m, n, pivsign;
33
34 /** Internal storage of pivot vector.
35 @serial pivot vector.
36 */
37 private int[] piv;
38
39 /* ------------------------
40 Constructor
41 * ------------------------ */
42
43 /** LU Decomposition
44 @param A Rectangular matrix
45 Structure to access L, U and piv.
46 */
47
48 public LUDecomposition(Matrix A) {
49
50 // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
51
52 LU = A.getArrayCopy();
53 m = A.getRowDimension();
54 n = A.getColumnDimension();
55 piv = new int[m];
56 for (int i = 0; i < m; i++) {
57 piv[i] = i;
58 }
59 pivsign = 1;
60 double[] LUrowi;
61 double[] LUcolj = new double[m];
62
63 // Outer loop.
64
65 for (int j = 0; j < n; j++) {
66
67 // Make a copy of the j-th column to localize references.
68
69 for (int i = 0; i < m; i++) {
70 LUcolj[i] = LU[i][j];
71 }
72
73 // Apply previous transformations.
74
75 for (int i = 0; i < m; i++) {
76 LUrowi = LU[i];
77
78 // Most of the time is spent in the following dot product.
79
80 int kmax = Math.min(i, j);
81 double s = 0.0;
82 for (int k = 0; k < kmax; k++) {
83 s += LUrowi[k] * LUcolj[k];
84 }
85
86 LUrowi[j] = LUcolj[i] -= s;
87 }
88
89 // Find pivot and exchange if necessary.
90
91 int p = j;
92 for (int i = j + 1; i < m; i++) {
93 if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
94 p = i;
95 }
96 }
97 if (p != j) {
98 for (int k = 0; k < n; k++) {
99 double t = LU[p][k];
100 LU[p][k] = LU[j][k];
101 LU[j][k] = t;
102 }
103 int k = piv[p];
104 piv[p] = piv[j];
105 piv[j] = k;
106 pivsign = -pivsign;
107 }
108
109 // Compute multipliers.
110
111 if (j < m & LU[j][j] != 0.0) {
112 for (int i = j + 1; i < m; i++) {
113 LU[i][j] /= LU[j][j];
114 }
115 }
116 }
117 }
118
119 /* ------------------------
120 Temporary, experimental code.
121 ------------------------ *\
122
123 \** LU Decomposition, computed by Gaussian elimination.
124 <P>
125 This constructor computes L and U with the "daxpy"-based elimination
126 algorithm used in LINPACK and MATLAB. In Java, we suspect the dot-product,
127 Crout algorithm will be faster. We have temporarily included this
128 constructor until timing experiments confirm this suspicion.
129 <P>
130 @param A Rectangular matrix
131 @param linpackflag Use Gaussian elimination. Actual value ignored.
132 @return Structure to access L, U and piv.
133 *\
134
135 public LUDecomposition (Matrix A, int linpackflag) {
136 // Initialize.
137 LU = A.getArrayCopy();
138 m = A.getRowDimension();
139 n = A.getColumnDimension();
140 piv = new int[m];
141 for (int i = 0; i < m; i++) {
142 piv[i] = i;
143 }
144 pivsign = 1;
145 // Main loop.
146 for (int k = 0; k < n; k++) {
147 // Find pivot.
148 int p = k;
149 for (int i = k+1; i < m; i++) {
150 if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) {
151 p = i;
152 }
153 }
154 // Exchange if necessary.
155 if (p != k) {
156 for (int j = 0; j < n; j++) {
157 double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t;
158 }
159 int t = piv[p]; piv[p] = piv[k]; piv[k] = t;
160 pivsign = -pivsign;
161 }
162 // Compute multipliers and eliminate k-th column.
163 if (LU[k][k] != 0.0) {
164 for (int i = k+1; i < m; i++) {
165 LU[i][k] /= LU[k][k];
166 for (int j = k+1; j < n; j++) {
167 LU[i][j] -= LU[i][k]*LU[k][j];
168 }
169 }
170 }
171 }
172 }
173
174 \* ------------------------
175 End of temporary code.
176 * ------------------------ */
177
178 /* ------------------------
179 Public Methods
180 * ------------------------ */
181
182 /** Is the matrix nonsingular?
183 @return true if U, and hence A, is nonsingular.
184 */
185
186 public boolean isNonsingular() {
187 for (int j = 0; j < n; j++) {
188 if (LU[j][j] == 0)
189 return false;
190 }
191 return true;
192 }
193
194 /** Return lower triangular factor
195 @return L
196 */
197
198 public Matrix getL() {
199 Matrix X = new Matrix(m, n);
200 double[][] L = X.getArray();
201 for (int i = 0; i < m; i++) {
202 for (int j = 0; j < n; j++) {
203 if (i > j) {
204 L[i][j] = LU[i][j];
205 } else if (i == j) {
206 L[i][j] = 1.0;
207 } else {
208 L[i][j] = 0.0;
209 }
210 }
211 }
212 return X;
213 }
214
215 /** Return upper triangular factor
216 @return U
217 */
218
219 public Matrix getU() {
220 Matrix X = new Matrix(n, n);
221 double[][] U = X.getArray();
222 for (int i = 0; i < n; i++) {
223 for (int j = 0; j < n; j++) {
224 if (i <= j) {
225 U[i][j] = LU[i][j];
226 } else {
227 U[i][j] = 0.0;
228 }
229 }
230 }
231 return X;
232 }
233
234 /** Return pivot permutation vector
235 @return piv
236 */
237
238 public int[] getPivot() {
239 int[] p = new int[m];
240 for (int i = 0; i < m; i++) {
241 p[i] = piv[i];
242 }
243 return p;
244 }
245
246 /** Return pivot permutation vector as a one-dimensional double array
247 @return (double) piv
248 */
249
250 public double[] getDoublePivot() {
251 double[] vals = new double[m];
252 for (int i = 0; i < m; i++) {
253 vals[i] = (double) piv[i];
254 }
255 return vals;
256 }
257
258 /** Determinant
259 @return det(A)
260 @exception IllegalArgumentException Matrix must be square
261 */
262
263 public double det() {
264 if (m != n) {
265 throw new IllegalArgumentException("Matrix must be square.");
266 }
267 double d = (double) pivsign;
268 for (int j = 0; j < n; j++) {
269 d *= LU[j][j];
270 }
271 return d;
272 }
273
274 /** Solve A*X = B
275 @param B A Matrix with as many rows as A and any number of columns.
276 @return X so that L*U*X = B(piv,:)
277 @exception IllegalArgumentException Matrix row dimensions must agree.
278 @exception RuntimeException Matrix is singular.
279 */
280
281 public Matrix solve(Matrix B) {
282 if (B.getRowDimension() != m) {
283 throw new IllegalArgumentException("Matrix row dimensions must agree.");
284 }
285 if (!this.isNonsingular()) {
286 throw new RuntimeException("Matrix is singular.");
287 }
288
289 // Copy right hand side with pivoting
290 int nx = B.getColumnDimension();
291 Matrix Xmat = B.getMatrix(piv, 0, nx - 1);
292 double[][] X = Xmat.getArray();
293
294 // Solve L*Y = B(piv,:)
295 for (int k = 0; k < n; k++) {
296 for (int i = k + 1; i < n; i++) {
297 for (int j = 0; j < nx; j++) {
298 X[i][j] -= X[k][j] * LU[i][k];
299 }
300 }
301 }
302 // Solve U*X = Y;
303 for (int k = n - 1; k >= 0; k--) {
304 for (int j = 0; j < nx; j++) {
305 X[k][j] /= LU[k][k];
306 }
307 for (int i = 0; i < k; i++) {
308 for (int j = 0; j < nx; j++) {
309 X[i][j] -= X[k][j] * LU[i][k];
310 }
311 }
312 }
313 return Xmat;
314 }
315 }
316