View Javadoc

1   package de.desy.acop.video.analysis;
2   
3   /**
4    * 
5    * <code>LUDecomposition</code> provides a method to solve a n-by-n system of linear
6    * equations. It makes a LU decomposition of the matrix (lower triangle L and upper 
7    * triangle U) and provides a solution to an arbitrary system A*x=b, where b
8    * is given by the user and x is returned. The implementation supports pivotting. 
9    *
10   * @author <a href="mailto:jaka.bobnar@cosylab.com">Jaka Bobnar</a>
11   *
12   */
13  public class LUDecomposition {
14  
15  	/** The matrix */
16  	private double[][] LU;
17  	/** Dimension */
18  	private int m;
19  	/** Pivot sign */
20  	private int pivsign;
21  	/** Pivot vector */
22  	private int[] pivot;
23  	 
24  
25  	/**
26  	 * Creates a new LUDecomposition object.
27  	 * 
28  	 * @param A the matrix to be decomposed
29  	 * @throws MathException if the matrix is singular
30  	 */
31  	public LUDecomposition(double[][] A) throws MathException {
32  		LU = getMatrixCopy(A);
33  		m = A.length;
34  		pivot = new int[m];
35  		for (int i = 0; i < m; i++) {
36  			pivot[i] = i;
37  		}
38  		pivsign = 1;
39  		decompose();
40  	}
41  
42  	/**
43  	 * Performs the LU decomposition of the matrix.
44  	 * 
45  	 * @throws MathException if the matrix is singular
46  	 */
47  	private void decompose() throws MathException {
48  		double[] LUrowi;
49  		double[] LUcolj = new double[m];
50  
51  		for (int j = 0; j < m; j++) {
52  			for (int i = 0; i < m; i++) {
53  				LUcolj[i] = LU[i][j];
54  			}
55  			for (int i = 0; i < m; i++) {
56  				LUrowi = LU[i];
57  				int kmax = Math.min(i,j);
58  				double s = 0.0;
59  				for (int k = 0; k < kmax; k++) {
60  					s += LUrowi[k] * LUcolj[k];
61  				}
62  				LUrowi[j] = LUcolj[i] -= s;
63  			}
64  			int p = j;
65  			for (int i = j + 1; i < m; i++) {
66  				if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
67  					p = i;
68  				}
69  			}
70  			if (p != j) {
71  				for (int k = 0; k < m; k++) {
72  					double t = LU[p][k];
73  					LU[p][k] = LU[j][k];
74  					LU[j][k] = t;
75  				}
76  				int k = pivot[p];
77  				pivot[p] = pivot[j];
78  				pivot[j] = k;
79  				pivsign = -pivsign;
80  			}
81  
82  			if (j < m & LU[j][j] != 0.0) {
83  				for (int i = j + 1; i < m; i++) {
84  					LU[i][j] /= LU[j][j];
85  				}
86  			}
87  		}
88  		if (isSingular()) {
89  			throw new MathException("Matrix is singular.");
90  		}
91  	}
92  
93  	/**
94  	 * Returns true if the given matrix is singular or false if not.
95  	 * 
96  	 * @return true if matrix is singular
97  	 */
98  	public boolean isSingular() {
99  		for (int j = 0; j < m; j++) {
100 			if (LU[j][j] == 0) return true;
101 		}
102 		return false;
103 	}
104 	
105 	/**
106 	 * Solves A*x = b
107 	 * 
108 	 * @param b vector of the dimension of this LU
109 	 * @return x so that L*U*x = b
110 	 * @throws MathException if the dimensions are not equal 
111 	 */
112 	public double[] solve(double[] b) throws MathException {
113 		if (b.length != m) {
114 			throw new MathException("Vector dimension must equal LU dimension (" + m +").");
115 		}
116 
117 		double[] x = new double[pivot.length];
118         for (int i = 0; i < x.length; i++) {
119        	 	x[i] = b[pivot[i]];
120         }
121 		
122 		// Solve L*Y = B(piv,:)
123 		for (int k = 0; k < m; k++) {
124 			for (int i = k + 1; i < m; i++) {
125 				x[i] -= x[k] * LU[i][k];
126 			}
127 		}
128 		// Solve U*X = Y;
129 		for (int k = m - 1; k >= 0; k--) {
130 			x[k] /= LU[k][k];
131 			for (int i = 0; i < k; i++) {
132 				x[i] -= x[k] * LU[i][k];
133 			}
134 		}
135 		return x;
136 	}
137 
138 	/**
139 	 * Returns the copy of the given matrix
140 	 * @param A the matrix to be copied 
141 	 * @return a copy of A
142 	 */
143 	private double[][] getMatrixCopy(double[][] A) {
144 		int m = A.length;
145 		double[][] C = new double[m][m];
146 		for (int i = 0; i < m; i++) {
147 			for (int j = 0; j < m; j++) {
148 				C[i][j] = A[i][j];
149 			}
150 		}
151 		return C;
152 	}
153 }