001/*
002 * To change this template, choose Tools | Templates
003 * and open the template in the editor.
004 */
005package org.jblas;
006
007import org.jblas.exceptions.LapackArgumentException;
008import org.jblas.exceptions.LapackPositivityException;
009import org.jblas.util.Permutations;
010import static org.jblas.util.Functions.min;
011
012/**
013 * Matrix which collects all kinds of decompositions.
014 */
015public 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}