public static SparseMatrix<Double> decomp(SparseMatrix<Double> a, boolean complete) { int m = a.rows(); SparseMatrix<Double> l = new SparseMatrix<Double>(m, m, Vectors.AS_SPARSE_REALS); //automatically initialzed to 0's for (int i = 0; i < m; i++) { if (complete) { for (int k = 0; k < (i + 1); k++) { chol(l, k, i, a); } } else { final Vector<Double> l_i = l.row(i); for (int k : l_i.keySet()) { if (k >= (i + 1)) { break; } chol(l, k, i, a); } } } return l; }
public static void chol(SparseMatrix<Double> l, int k, int i, SparseMatrix<Double> a) throws IllegalArgumentException { double sum = 0; final Vector<Double> l_k = l.row(k); for (int j : l_k.keySet()) { sum += l.doubleValue(i, j) * l_k.doubleValue(j); } double a_ii = a.doubleValue(i, i); if (i == k) { if (a_ii - sum < 0) { throw new IllegalArgumentException("Matrix not positive definite"); } l.set(i, k, Math.sqrt(a_ii - sum)); } else { l.set(i, k, (a.doubleValue(i, k) - sum) / l.doubleValue(k, k)); } }
@Override public <M extends Number> void sub(Matrix<M> matrix) { assert (m == matrix.rows()); assert (n == matrix.cols()); if (matrix instanceof SparseMatrix) { final Vector<M>[] arr2 = ((SparseMatrix<M>) matrix).arr; for (int i = 0; i < m; i++) { if (arr[i] != null && arr2 != null) { arr[i].sub(arr2[i]); } else if (arr2 != null) { arr[i] = using.make(n, defaultValue); arr[i].sub(arr2[i]); } } } else { for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { row(i).sub(j, matrix.doubleValue(i, j)); } } } }
@Override public <M extends Number> void add(Matrix<M> matrix) { assert (m == matrix.rows()); assert (n == matrix.cols()); if (matrix instanceof SparseMatrix) { final Vector<M>[] arr2 = ((SparseMatrix<M>) matrix).arr; for (int i = 0; i < m; i++) { if (arr[i] != null && arr2 != null) { arr[i].add(arr2[i]); } else if (arr2 != null) { arr[i] = using.make(n, defaultValue); arr[i].add(arr2[i]); } } } else { for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { row(i).add(j, matrix.doubleValue(i, j)); } } } }