001    /*
002     * To change this template, choose Tools | Templates
003     * and open the template in the editor.
004     */
005    package org.jblas;
006    
007    import org.jblas.exceptions.LapackArgumentException;
008    import org.jblas.exceptions.LapackPositivityException;
009    import org.jblas.util.Permutations;
010    import static org.jblas.util.Functions.min;
011    
012    /**
013     * Matrix which collects all kinds of decompositions.
014     */
015    public class Decompose {
016    
017        public static class LUDecomposition<T> {
018    
019            public T l;
020            public T u;
021            public T p;
022    
023            public LUDecomposition(T l, T u, T p) {
024                this.l = l;
025                this.u = u;
026                this.p = p;
027            }
028        }
029    
030        public static LUDecomposition<DoubleMatrix> lu(DoubleMatrix A) {
031            int[] ipiv = new int[min(A.rows, A.columns)];
032            DoubleMatrix result = A.dup();
033            NativeBlas.dgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0);
034    
035            // collect result
036            DoubleMatrix l = new DoubleMatrix(A.rows, min(A.rows, A.columns));
037            DoubleMatrix u = new DoubleMatrix(min(A.columns, A.rows), A.columns);
038            decomposeLowerUpper(result, l, u);
039            DoubleMatrix p = Permutations.permutationMatrixFromPivotIndices(A.rows, ipiv);
040            return new LUDecomposition<DoubleMatrix>(l, u, p);
041        }
042    
043        private static void decomposeLowerUpper(DoubleMatrix A, DoubleMatrix L, DoubleMatrix U) {
044            for (int i = 0; i < A.rows; i++) {
045                for (int j = 0; j < A.columns; j++) {
046                    if (i < j) {
047                        U.put(i, j, A.get(i, j));
048                    } else if (i == j) {
049                        U.put(i, i, A.get(i, i));
050                        L.put(i, i, 1.0);
051                    } else {
052                        L.put(i, j, A.get(i, j));
053                    }
054    
055                }
056            }
057        }
058    
059        /**
060         * Compute Cholesky decomposition of A
061         *
062         * @param A symmetric, positive definite matrix (only upper half is used)
063         * @return upper triangular matrix U such that  A = U' * U
064         */
065        public static DoubleMatrix cholesky(DoubleMatrix A) {
066            DoubleMatrix result = A.dup();
067            int info = NativeBlas.dpotrf('U', A.rows, result.data, 0, A.rows);
068            if (info < 0) {
069                throw new LapackArgumentException("DPOTRF", -info);
070            } else if (info > 0) {
071                throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite.");
072            }
073            clearLower(result);
074            return result;
075        }
076    
077        private static void clearLower(DoubleMatrix A) {
078            for (int j = 0; j < A.columns; j++)
079                for (int i = j + 1; i < A.rows; i++)
080                    A.put(i, j, 0.0);
081        }
082    }