1 package Jama;
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 @version $Id: LUDecomposition.java,v 1.1.1.1 2010/11/30 21:31:59 jeremy Exp $
15 */
16
17 public class LUDecomposition implements java.io.Serializable {
18
19 /* ------------------------
20 Class variables
21 * ------------------------ */
22
23 /** Array for internal storage of decomposition.
24 @serial internal array storage.
25 */
26 private double[][] LU;
27
28 /** Row and column dimensions, and pivot sign.
29 @serial column dimension.
30 @serial row dimension.
31 @serial pivot sign.
32 */
33 private int m, n, pivsign;
34
35 /** Internal storage of pivot vector.
36 @serial pivot vector.
37 */
38 private int[] piv;
39
40 /* ------------------------
41 Constructor
42 * ------------------------ */
43
44 /** LU Decomposition
45 @param A Rectangular matrix
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]; LU[p][k] = LU[j][k]; LU[j][k] = t;
100 }
101 int k = piv[p]; piv[p] = piv[j]; piv[j] = k;
102 pivsign = -pivsign;
103 }
104
105 // Compute multipliers.
106
107 if (j < m & LU[j][j] != 0.0) {
108 for (int i = j+1; i < m; i++) {
109 LU[i][j] /= LU[j][j];
110 }
111 }
112 }
113 }
114
115 /* ------------------------
116 Temporary, experimental code.
117 ------------------------ *\
118
119 \** LU Decomposition, computed by Gaussian elimination.
120 <P>
121 This constructor computes L and U with the "daxpy"-based elimination
122 algorithm used in LINPACK and MATLAB. In Java, we suspect the dot-product,
123 Crout algorithm will be faster. We have temporarily included this
124 constructor until timing experiments confirm this suspicion.
125 <P>
126 @param A Rectangular matrix
127 @param linpackflag Use Gaussian elimination. Actual value ignored.
128 @return Structure to access L, U and piv.
129 *\
130
131 public LUDecomposition (Matrix A, int linpackflag) {
132 // Initialize.
133 LU = A.getArrayCopy();
134 m = A.getRowDimension();
135 n = A.getColumnDimension();
136 piv = new int[m];
137 for (int i = 0; i < m; i++) {
138 piv[i] = i;
139 }
140 pivsign = 1;
141 // Main loop.
142 for (int k = 0; k < n; k++) {
143 // Find pivot.
144 int p = k;
145 for (int i = k+1; i < m; i++) {
146 if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) {
147 p = i;
148 }
149 }
150 // Exchange if necessary.
151 if (p != k) {
152 for (int j = 0; j < n; j++) {
153 double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t;
154 }
155 int t = piv[p]; piv[p] = piv[k]; piv[k] = t;
156 pivsign = -pivsign;
157 }
158 // Compute multipliers and eliminate k-th column.
159 if (LU[k][k] != 0.0) {
160 for (int i = k+1; i < m; i++) {
161 LU[i][k] /= LU[k][k];
162 for (int j = k+1; j < n; j++) {
163 LU[i][j] -= LU[i][k]*LU[k][j];
164 }
165 }
166 }
167 }
168 }
169
170 \* ------------------------
171 End of temporary code.
172 * ------------------------ */
173
174 /* ------------------------
175 Public Methods
176 * ------------------------ */
177
178 /** Is the matrix nonsingular?
179 @return true if U, and hence A, is nonsingular.
180 */
181
182 public boolean isNonsingular () {
183 for (int j = 0; j < n; j++) {
184 if (LU[j][j] == 0)
185 return false;
186 }
187 return true;
188 }
189
190 /** Return lower triangular factor
191 @return L
192 */
193
194 public Matrix getL () {
195 Matrix X = new Matrix(m,n);
196 double[][] L = X.getArray();
197 for (int i = 0; i < m; i++) {
198 for (int j = 0; j < n; j++) {
199 if (i > j) {
200 L[i][j] = LU[i][j];
201 } else if (i == j) {
202 L[i][j] = 1.0;
203 } else {
204 L[i][j] = 0.0;
205 }
206 }
207 }
208 return X;
209 }
210
211 /** Return upper triangular factor
212 @return U
213 */
214
215 public Matrix getU () {
216 Matrix X = new Matrix(n,n);
217 double[][] U = X.getArray();
218 for (int i = 0; i < n; i++) {
219 for (int j = 0; j < n; j++) {
220 if (i <= j) {
221 U[i][j] = LU[i][j];
222 } else {
223 U[i][j] = 0.0;
224 }
225 }
226 }
227 return X;
228 }
229
230 /** Return pivot permutation vector
231 @return piv
232 */
233
234 public int[] getPivot () {
235 int[] p = new int[m];
236 for (int i = 0; i < m; i++) {
237 p[i] = piv[i];
238 }
239 return p;
240 }
241
242 /** Return pivot permutation vector as a one-dimensional double array
243 @return (double) piv
244 */
245
246 public double[] getDoublePivot () {
247 double[] vals = new double[m];
248 for (int i = 0; i < m; i++) {
249 vals[i] = (double) piv[i];
250 }
251 return vals;
252 }
253
254 /** Determinant
255 @return det(A)
256 @exception IllegalArgumentException Matrix must be square
257 */
258
259 public double det () {
260 if (m != n) {
261 throw new IllegalArgumentException("Matrix must be square.");
262 }
263 double d = (double) pivsign;
264 for (int j = 0; j < n; j++) {
265 d *= LU[j][j];
266 }
267 return d;
268 }
269
270 /** Solve A*X = B
271 @param B A Matrix with as many rows as A and any number of columns.
272 @return X so that L*U*X = B(piv,:)
273 @exception IllegalArgumentException Matrix row dimensions must agree.
274 @exception RuntimeException Matrix is singular.
275 */
276
277 public Matrix solve (Matrix B) {
278 if (B.getRowDimension() != m) {
279 throw new IllegalArgumentException("Matrix row dimensions must agree.");
280 }
281 if (!this.isNonsingular()) {
282 throw new RuntimeException("Matrix is singular.");
283 }
284
285 // Copy right hand side with pivoting
286 int nx = B.getColumnDimension();
287 Matrix Xmat = B.getMatrix(piv,0,nx-1);
288 double[][] X = Xmat.getArray();
289
290 // Solve L*Y = B(piv,:)
291 for (int k = 0; k < n; k++) {
292 for (int i = k+1; i < n; i++) {
293 for (int j = 0; j < nx; j++) {
294 X[i][j] -= X[k][j]*LU[i][k];
295 }
296 }
297 }
298 // Solve U*X = Y;
299 for (int k = n-1; k >= 0; k--) {
300 for (int j = 0; j < nx; j++) {
301 X[k][j] /= LU[k][k];
302 }
303 for (int i = 0; i < k; i++) {
304 for (int j = 0; j < nx; j++) {
305 X[i][j] -= X[k][j]*LU[i][k];
306 }
307 }
308 }
309 return Xmat;
310 }
311 }