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     * Class to hold an LU decomposition result.
018     *
019     * Contains a lower matrix L, and upper matrix U, and a permutation matrix
020     * P such that P*L*U is the original matrix.
021     * @param <T>
022     */
023    public static class LUDecomposition<T> {
024
025        public T l;
026        public T u;
027        public T p;
028
029        public LUDecomposition(T l, T u, T p) {
030            this.l = l;
031            this.u = u;
032            this.p = p;
033        }
034
035        @Override
036        public String toString() {
037          return String.format("<LUDecomposition L=%s U=%s P=%s>", l, u, p);
038        }
039    }
040
041    /**
042     * Compute LU Decomposition of a general matrix.
043     *
044     * Computes the LU decomposition using GETRF. Returns three matrices L, U, P,
045     * where L is lower diagonal, U is upper diagonal, and P is a permutation
046     * matrix such that A = P * L * U.
047     *
048     * @param A general matrix
049     * @return An LUDecomposition object.
050     */
051    public static LUDecomposition<DoubleMatrix> lu(DoubleMatrix A) {
052        int[] ipiv = new int[min(A.rows, A.columns)];
053        DoubleMatrix result = A.dup();
054        NativeBlas.dgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0);
055
056        // collect result
057        DoubleMatrix l = new DoubleMatrix(A.rows, min(A.rows, A.columns));
058        DoubleMatrix u = new DoubleMatrix(min(A.columns, A.rows), A.columns);
059        decomposeLowerUpper(result, l, u);
060        DoubleMatrix p = Permutations.permutationDoubleMatrixFromPivotIndices(A.rows, ipiv);
061        return new LUDecomposition<DoubleMatrix>(l, u, p);
062    }
063
064    private static void decomposeLowerUpper(DoubleMatrix A, DoubleMatrix L, DoubleMatrix U) {
065        for (int i = 0; i < A.rows; i++) {
066            for (int j = 0; j < A.columns; j++) {
067                if (i < j) {
068                    U.put(i, j, A.get(i, j));
069                } else if (i == j) {
070                    U.put(i, i, A.get(i, i));
071                    L.put(i, i, 1.0);
072                } else {
073                    L.put(i, j, A.get(i, j));
074                }
075
076            }
077        }
078    }
079
080    /**if (info )
081     * Compute Cholesky decomposition of A
082     *
083     * @param A symmetric, positive definite matrix (only upper half is used)
084     * @return upper triangular matrix U such that  A = U' * U
085     */
086    public static FloatMatrix cholesky(FloatMatrix A) {
087        FloatMatrix result = A.dup();
088        int info = NativeBlas.spotrf('U', A.rows, result.data, 0, A.rows);
089        if (info < 0) {
090            throw new LapackArgumentException("DPOTRF", -info);
091        } else if (info > 0) {
092            throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite.");
093        }
094        clearLower(result);
095        return result;
096    }
097
098    private static void clearLower(FloatMatrix A) {
099        for (int j = 0; j < A.columns; j++)
100            for (int i = j + 1; i < A.rows; i++)
101                A.put(i, j, 0.0f);
102    }
103
104  /**
105   * Compute LU Decomposition of a general matrix.
106   *
107   * Computes the LU decomposition using GETRF. Returns three matrices L, U, P,
108   * where L is lower diagonal, U is upper diagonal, and P is a permutation
109   * matrix such that A = P * L * U.
110   *
111   * @param A general matrix
112   * @return An LUDecomposition object.
113   */
114  public static LUDecomposition<FloatMatrix> lu(FloatMatrix A) {
115      int[] ipiv = new int[min(A.rows, A.columns)];
116      FloatMatrix result = A.dup();
117      NativeBlas.sgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0);
118
119      // collect result
120      FloatMatrix l = new FloatMatrix(A.rows, min(A.rows, A.columns));
121      FloatMatrix u = new FloatMatrix(min(A.columns, A.rows), A.columns);
122      decomposeLowerUpper(result, l, u);
123      FloatMatrix p = Permutations.permutationFloatMatrixFromPivotIndices(A.rows, ipiv);
124      return new LUDecomposition<FloatMatrix>(l, u, p);
125  }
126
127  private static void decomposeLowerUpper(FloatMatrix A, FloatMatrix L, FloatMatrix U) {
128      for (int i = 0; i < A.rows; i++) {
129          for (int j = 0; j < A.columns; j++) {
130              if (i < j) {
131                  U.put(i, j, A.get(i, j));
132              } else if (i == j) {
133                  U.put(i, i, A.get(i, i));
134                  L.put(i, i, 1.0f);
135              } else {
136                  L.put(i, j, A.get(i, j));
137              }
138
139          }
140      }
141  }
142
143  /**
144   * Compute Cholesky decomposition of A
145   *
146   * @param A symmetric, positive definite matrix (only upper half is used)
147   * @return upper triangular matrix U such that  A = U' * U
148   */
149  public static DoubleMatrix cholesky(DoubleMatrix A) {
150      DoubleMatrix result = A.dup();
151      int info = NativeBlas.dpotrf('U', A.rows, result.data, 0, A.rows);
152      if (info < 0) {
153          throw new LapackArgumentException("DPOTRF", -info);
154      } else if (info > 0) {
155          throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite.");
156      }
157      clearLower(result);
158      return result;
159  }
160
161  private static void clearLower(DoubleMatrix A) {
162      for (int j = 0; j < A.columns; j++)
163          for (int i = j + 1; i < A.rows; i++)
164              A.put(i, j, 0.0);
165  }
166
167  /**
168   * Class to represent a QR decomposition.
169   *
170   * @param <T>
171   */
172  public static class QRDecomposition<T> {
173    public T q;
174    public T r;
175
176    QRDecomposition(T q, T r) {
177      this.q = q;
178      this.r = r;
179    }
180
181    @Override
182    public String toString() {
183      return "<Q=" + q + " R=" + r + ">";
184    }
185  }
186
187  /**
188   * QR decomposition.
189   *
190   * Decomposes (m,n) matrix A into a (m,m) matrix Q and an (m,n) matrix R such that
191   * Q is orthogonal, R is upper triangular and Q * R = A
192   *
193   * Note that if A has more rows than columns, then the lower rows of R will contain
194   * only zeros, such that the corresponding later columns of Q do not enter the computation
195   * at all. For some reason, LAPACK does not properly normalize those columns.
196   *
197   * @param A matrix
198   * @return QR decomposition
199   */
200  public static QRDecomposition<DoubleMatrix> qr(DoubleMatrix A) {
201    int minmn = min(A.rows, A.columns);
202    DoubleMatrix result = A.dup();
203    DoubleMatrix tau = new DoubleMatrix(minmn);
204    SimpleBlas.geqrf(result, tau);
205    DoubleMatrix R = new DoubleMatrix(A.rows, A.columns);
206    for (int i = 0; i < A.rows; i++) {
207      for (int j = i; j < A.columns; j++) {
208        R.put(i, j, result.get(i, j));
209      }
210    }
211    DoubleMatrix Q = DoubleMatrix.eye(A.rows);
212    SimpleBlas.ormqr('L', 'N', result, tau, Q);
213    return new QRDecomposition<DoubleMatrix>(Q, R);
214  }
215  
216  /**
217   * QR decomposition.
218   *
219   * Decomposes (m,n) matrix A into a (m,m) matrix Q and an (m,n) matrix R such that
220   * Q is orthogonal, R is upper triangular and Q * R = A
221   *
222   * @param A matrix
223   * @return QR decomposition
224   */
225  public static QRDecomposition<FloatMatrix> qr(FloatMatrix A) {
226    int minmn = min(A.rows, A.columns);
227    FloatMatrix result = A.dup();
228    FloatMatrix tau = new FloatMatrix(minmn);
229    SimpleBlas.geqrf(result, tau);
230    FloatMatrix R = new FloatMatrix(A.rows, A.columns);
231    for (int i = 0; i < A.rows; i++) {
232      for (int j = i; j < A.columns; j++) {
233        R.put(i, j, result.get(i, j));
234      }
235    }
236    FloatMatrix Q = FloatMatrix.eye(A.rows);
237    SimpleBlas.ormqr('L', 'N', result, tau, Q);
238    return new QRDecomposition<FloatMatrix>(Q, R);
239  }
240}