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