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 public static class LUDecomposition<T> { 018 019 public T l; 020 public T u; 021 public T p; 022 023 public LUDecomposition(T l, T u, T p) { 024 this.l = l; 025 this.u = u; 026 this.p = p; 027 } 028 } 029 030 public static LUDecomposition<DoubleMatrix> lu(DoubleMatrix A) { 031 int[] ipiv = new int[min(A.rows, A.columns)]; 032 DoubleMatrix result = A.dup(); 033 NativeBlas.dgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0); 034 035 // collect result 036 DoubleMatrix l = new DoubleMatrix(A.rows, min(A.rows, A.columns)); 037 DoubleMatrix u = new DoubleMatrix(min(A.columns, A.rows), A.columns); 038 decomposeLowerUpper(result, l, u); 039 DoubleMatrix p = Permutations.permutationMatrixFromPivotIndices(A.rows, ipiv); 040 return new LUDecomposition<DoubleMatrix>(l, u, p); 041 } 042 043 private static void decomposeLowerUpper(DoubleMatrix A, DoubleMatrix L, DoubleMatrix U) { 044 for (int i = 0; i < A.rows; i++) { 045 for (int j = 0; j < A.columns; j++) { 046 if (i < j) { 047 U.put(i, j, A.get(i, j)); 048 } else if (i == j) { 049 U.put(i, i, A.get(i, i)); 050 L.put(i, i, 1.0); 051 } else { 052 L.put(i, j, A.get(i, j)); 053 } 054 055 } 056 } 057 } 058 059 /** 060 * Compute Cholesky decomposition of A 061 * 062 * @param A symmetric, positive definite matrix (only upper half is used) 063 * @return upper triangular matrix U such that A = U' * U 064 */ 065 public static DoubleMatrix cholesky(DoubleMatrix A) { 066 DoubleMatrix result = A.dup(); 067 int info = NativeBlas.dpotrf('U', A.rows, result.data, 0, A.rows); 068 if (info < 0) { 069 throw new LapackArgumentException("DPOTRF", -info); 070 } else if (info > 0) { 071 throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite."); 072 } 073 clearLower(result); 074 return result; 075 } 076 077 private static void clearLower(DoubleMatrix A) { 078 for (int j = 0; j < A.columns; j++) 079 for (int i = j + 1; i < A.rows; i++) 080 A.put(i, j, 0.0); 081 } 082 }