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)); } }
public static void denseInplaceDecomp(SparseMatrix<Double> a) { int n = a.rows(); int i, j, k; double sum; double[] diagonal = new double[n]; for (i = 0; i < n; i++) { System.err.print("."); for (j = i; j < n; j++) { for (sum = a.doubleValue(i, j), k = i - 1; k >= 0; k--) { sum -= a.doubleValue(i, k) * a.doubleValue(j, k); } if (i == j) { if (sum <= 0.0) { throw new IllegalArgumentException("Matrix not positive definite"); } diagonal[i] = Math.sqrt(sum); } else { a.set(j, i, sum / diagonal[i]); a.set(i, j, 0); } } } for (i = 0; i < n; i++) { a.set(i, i, diagonal[i]); } }
public static Vector<Double> solve(SparseMatrix<Double> a, Vector<Double> b) { assert (a.cols() == b.length()); final int N = b.length(); Vector<Double> y = new SparseRealArray(N); for (int i = 0; i < N; i++) { double sum = b.doubleValue(i); for (int j = 0; j < i; j++) { sum -= a.doubleValue(i, j) * y.doubleValue(j); } y.put(i, sum / a.doubleValue(i, i)); } return y; }
public static Vector<Double> solveT(SparseMatrix<Double> a, Vector<Double> b) { assert (a.cols() == b.length()); final int N = b.length(); Vector<Double> y = new SparseRealArray(N); for (int i = N - 1; i >= 0; i--) { double sum = b.doubleValue(i); for (int j = i + 1; j < N; j++) { sum -= a.doubleValue(j, i) * y.doubleValue(j); } y.put(i, sum / a.doubleValue(i, i)); } return y; }
@Override public double[][] toDoubleArray() { double[][] d = new double[rows()][cols()]; for (int i = 0; i < rows(); i++) { for (int j = 0; j < rows(); j++) { d[i][j] = doubleValue(i, j); } } return d; } }