Source code

001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.structure.jama;
022
023        /** LU Decomposition.
024        <P>
025        For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
026        unit lower triangular matrix L, an n-by-n upper triangular matrix U,
027        and a permutation vector piv of length m so that A(piv,:) = L*U.
028        If m < n, then L is m-by-m and U is m-by-n.
029        <P>
030        The LU decompostion with pivoting always exists, even if the matrix is
031        singular, so the constructor will never fail.  The primary use of the
032        LU decomposition is in the solution of square systems of simultaneous
033        linear equations.  This will fail if isNonsingular() returns false.
034        */
035
036public class LUDecomposition implements java.io.Serializable {
037
038         static final long serialVersionUID = 9271028462937843l;
039
040/* ------------------------
041        Class variables
042 * ------------------------ */
043
044        /** Array for internal storage of decomposition.
045        @serial internal array storage.
046        */
047        private double[][] LU;
048
049        /** Row and column dimensions, and pivot sign.
050        @serial column dimension.
051        @serial row dimension.
052        @serial pivot sign.
053        */
054        private int m, n, pivsign;
055
056        /** Internal storage of pivot vector.
057        @serial pivot vector.
058        */
059        private int[] piv;
060
061/* ------------------------
062        Constructor
063 * ------------------------ */
064
065        /** LU Decomposition provides a data structure  to access L, U and piv.
066        @param  A   Rectangular matrix
067        */
068
069        public LUDecomposition (Matrix A) {
070
071        // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
072
073                LU = A.getArrayCopy();
074                m = A.getRowDimension();
075                n = A.getColumnDimension();
076                piv = new int[m];
077                for (int i = 0; i < m; i++) {
078                        piv[i] = i;
079                }
080                pivsign = 1;
081                double[] LUrowi;
082                double[] LUcolj = new double[m];
083
084                // Outer loop.
085
086                for (int j = 0; j < n; j++) {
087
088                        // Make a copy of the j-th column to localize references.
089
090                        for (int i = 0; i < m; i++) {
091                                LUcolj[i] = LU[i][j];
092                        }
093
094                        // Apply previous transformations.
095
096                        for (int i = 0; i < m; i++) {
097                                LUrowi = LU[i];
098
099                                // Most of the time is spent in the following dot product.
100
101                                int kmax = Math.min(i,j);
102                                double s = 0.0;
103                                for (int k = 0; k < kmax; k++) {
104                                        s += LUrowi[k]*LUcolj[k];
105                                }
106
107                                LUrowi[j] = LUcolj[i] -= s;
108                        }
109
110                        // Find pivot and exchange if necessary.
111
112                        int p = j;
113                        for (int i = j+1; i < m; i++) {
114                                if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
115                                        p = i;
116                                }
117                        }
118                        if (p != j) {
119                                for (int k = 0; k < n; k++) {
120                                        double t = LU[p][k]; LU[p][k] = LU[j][k]; LU[j][k] = t;
121                                }
122                                int k = piv[p]; piv[p] = piv[j]; piv[j] = k;
123                                pivsign = -pivsign;
124                        }
125
126                        // Compute multipliers.
127
128                        if (j < m && LU[j][j] != 0.0) {
129                                for (int i = j+1; i < m; i++) {
130                                        LU[i][j] /= LU[j][j];
131                                }
132                        }
133                }
134        }
135
136/* ------------------------
137        Temporary, experimental code.
138        ------------------------ *\
139
140        \** LU Decomposition, computed by Gaussian elimination.
141        <P>
142        This constructor computes L and U with the "daxpy"-based elimination
143        algorithm used in LINPACK and MATLAB.  In Java, we suspect the dot-product,
144        Crout algorithm will be faster.  We have temporarily included this
145        constructor until timing experiments confirm this suspicion.
146        <P>
147        @param  A             Rectangular matrix
148        @param  linpackflag   Use Gaussian elimination.  Actual value ignored.
149        @return               Structure to access L, U and piv.
150        *\
151
152        public LUDecomposition (Matrix A, int linpackflag) {
153                // Initialize.
154                LU = A.getArrayCopy();
155                m = A.getRowDimension();
156                n = A.getColumnDimension();
157                piv = new int[m];
158                for (int i = 0; i < m; i++) {
159                        piv[i] = i;
160                }
161                pivsign = 1;
162                // Main loop.
163                for (int k = 0; k < n; k++) {
164                        // Find pivot.
165                        int p = k;
166                        for (int i = k+1; i < m; i++) {
167                                if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) {
168                                        p = i;
169                                }
170                        }
171                        // Exchange if necessary.
172                        if (p != k) {
173                                for (int j = 0; j < n; j++) {
174                                        double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t;
175                                }
176                                int t = piv[p]; piv[p] = piv[k]; piv[k] = t;
177                                pivsign = -pivsign;
178                        }
179                        // Compute multipliers and eliminate k-th column.
180                        if (LU[k][k] != 0.0) {
181                                for (int i = k+1; i < m; i++) {
182                                        LU[i][k] /= LU[k][k];
183                                        for (int j = k+1; j < n; j++) {
184                                                LU[i][j] -= LU[i][k]*LU[k][j];
185                                        }
186                                }
187                        }
188                }
189        }
190
191\* ------------------------
192        End of temporary code.
193 * ------------------------ */
194
195/* ------------------------
196        Public Methods
197 * ------------------------ */
198
199        /** Is the matrix nonsingular?
200        @return     true if U, and hence A, is nonsingular.
201        */
202
203        public boolean isNonsingular () {
204                for (int j = 0; j < n; j++) {
205                        if (LU[j][j] == 0)
206                                return false;
207                }
208                return true;
209        }
210
211        /** Return lower triangular factor
212        @return     L
213        */
214
215        public Matrix getL () {
216                Matrix X = new Matrix(m,n);
217                double[][] L = X.getArray();
218                for (int i = 0; i < m; i++) {
219                        for (int j = 0; j < n; j++) {
220                                if (i > j) {
221                                        L[i][j] = LU[i][j];
222                                } else if (i == j) {
223                                        L[i][j] = 1.0;
224                                } else {
225                                        L[i][j] = 0.0;
226                                }
227                        }
228                }
229                return X;
230        }
231
232        /** Return upper triangular factor
233        @return     U
234        */
235
236        public Matrix getU () {
237                Matrix X = new Matrix(n,n);
238                double[][] U = X.getArray();
239                for (int i = 0; i < n; i++) {
240                        for (int j = 0; j < n; j++) {
241                                if (i <= j) {
242                                        U[i][j] = LU[i][j];
243                                } else {
244                                        U[i][j] = 0.0;
245                                }
246                        }
247                }
248                return X;
249        }
250
251        /** Return pivot permutation vector
252        @return     piv
253        */
254
255        public int[] getPivot () {
256                int[] p = new int[m];
257                for (int i = 0; i < m; i++) {
258                        p[i] = piv[i];
259                }
260                return p;
261        }
262
263        /** Return pivot permutation vector as a one-dimensional double array
264        @return     (double) piv
265        */
266
267        public double[] getDoublePivot () {
268                double[] vals = new double[m];
269                for (int i = 0; i < m; i++) {
270                        vals[i] = piv[i];
271                }
272                return vals;
273        }
274
275        /** Determinant
276        @return     det(A)
277        @exception  IllegalArgumentException  Matrix must be square
278        */
279
280        public double det () {
281                if (m != n) {
282                        throw new IllegalArgumentException("Matrix must be square.");
283                }
284                double d = pivsign;
285                for (int j = 0; j < n; j++) {
286                        d *= LU[j][j];
287                }
288                return d;
289        }
290
291        /** Solve A*X = B
292        @param  B   A Matrix with as many rows as A and any number of columns.
293        @return     X so that L*U*X = B(piv,:)
294        @exception  IllegalArgumentException Matrix row dimensions must agree.
295        @exception  RuntimeException  Matrix is singular.
296        */
297
298        public Matrix solve (Matrix B) {
299                if (B.getRowDimension() != m) {
300                        throw new IllegalArgumentException("Matrix row dimensions must agree.");
301                }
302                if (!this.isNonsingular()) {
303                        throw new RuntimeException("Matrix is singular.");
304                }
305
306                // Copy right hand side with pivoting
307                int nx = B.getColumnDimension();
308                Matrix Xmat = B.getMatrix(piv,0,nx-1);
309                double[][] X = Xmat.getArray();
310
311                // Solve L*Y = B(piv,:)
312                for (int k = 0; k < n; k++) {
313                        for (int i = k+1; i < n; i++) {
314                                for (int j = 0; j < nx; j++) {
315                                        X[i][j] -= X[k][j]*LU[i][k];
316                                }
317                        }
318                }
319                // Solve U*X = Y;
320                for (int k = n-1; k >= 0; k--) {
321                        for (int j = 0; j < nx; j++) {
322                                X[k][j] /= LU[k][k];
323                        }
324                        for (int i = 0; i < k; i++) {
325                                for (int j = 0; j < nx; j++) {
326                                        X[i][j] -= X[k][j]*LU[i][k];
327                                }
328                        }
329                }
330                return Xmat;
331        }
332}