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 }