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        /**
018         * Class to hold an LU decomposition result.
019         *
020         * Contains a lower matrix L, and upper matrix U, and a permutation matrix
021         * P such that P*L*U is the original matrix.
022         * @param <T>
023         */
024        public static class LUDecomposition<T> {
025    
026            public T l;
027            public T u;
028            public T p;
029    
030            public LUDecomposition(T l, T u, T p) {
031                this.l = l;
032                this.u = u;
033                this.p = p;
034            }
035        }
036    
037        /**
038         * Compute LU Decomposition of a general matrix.
039         *
040         * Computes the LU decomposition using GETRF. Returns three matrices L, U, P,
041         * where L is lower diagonal, U is upper diagonal, and P is a permutation
042         * matrix such that A = P * L * U.
043         *
044         * @param A general matrix
045         * @return An LUDecomposition object.
046         */
047        public static LUDecomposition<DoubleMatrix> lu(DoubleMatrix A) {
048            int[] ipiv = new int[min(A.rows, A.columns)];
049            DoubleMatrix result = A.dup();
050            NativeBlas.dgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0);
051    
052            // collect result
053            DoubleMatrix l = new DoubleMatrix(A.rows, min(A.rows, A.columns));
054            DoubleMatrix u = new DoubleMatrix(min(A.columns, A.rows), A.columns);
055            decomposeLowerUpper(result, l, u);
056            DoubleMatrix p = Permutations.permutationMatrixFromPivotIndices(A.rows, ipiv);
057            return new LUDecomposition<DoubleMatrix>(l, u, p);
058        }
059    
060        private static void decomposeLowerUpper(DoubleMatrix A, DoubleMatrix L, DoubleMatrix U) {
061            for (int i = 0; i < A.rows; i++) {
062                for (int j = 0; j < A.columns; j++) {
063                    if (i < j) {
064                        U.put(i, j, A.get(i, j));
065                    } else if (i == j) {
066                        U.put(i, i, A.get(i, i));
067                        L.put(i, i, 1.0);
068                    } else {
069                        L.put(i, j, A.get(i, j));
070                    }
071    
072                }
073            }
074        }
075    
076        /**
077         * Compute Cholesky decomposition of A
078         *
079         * @param A symmetric, positive definite matrix (only upper half is used)
080         * @return upper triangular matrix U such that  A = U' * U
081         */
082        public static DoubleMatrix cholesky(DoubleMatrix A) {
083            DoubleMatrix result = A.dup();
084            int info = NativeBlas.dpotrf('U', A.rows, result.data, 0, A.rows);
085            if (info < 0) {
086                throw new LapackArgumentException("DPOTRF", -info);
087            } else if (info > 0) {
088                throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite.");
089            }
090            clearLower(result);
091            return result;
092        }
093    
094        private static void clearLower(DoubleMatrix A) {
095            for (int j = 0; j < A.columns; j++)
096                for (int i = j + 1; i < A.rows; i++)
097                    A.put(i, j, 0.0);
098        }
099    }