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