private DoubleMatrix1D upperSolve(DoubleMatrix1D b, DoubleMatrix1D x) { double[] bd = ((DenseDoubleMatrix1D) b).elements(); double[] xd = ((DenseDoubleMatrix1D) x).elements(); int rows = LU.rows(); for (int i = rows - 1; i >= 0; --i) { // Get row i SparseDoubleMatrix1D row = LU.viewRow(i); int used = (int) row.size(); // xi = (bi - sum[j>i] Uij * xj) / Uii double sum = 0; for (int j = i + 1; j < used; ++j) sum += row.getQuick(j) * xd[j]; xd[i] = (bd[i] - sum) / row.getQuick(i); } return x; }
private DoubleMatrix1D upperSolve(DoubleMatrix1D b, DoubleMatrix1D x) { double[] bd = ((DenseDoubleMatrix1D) b).elements(); double[] xd = ((DenseDoubleMatrix1D) x).elements(); int rows = LU.rows(); for (int i = rows - 1; i >= 0; --i) { // Get row i SparseDoubleMatrix1D row = LU.viewRow(i); int used = (int) row.size(); // xi = (bi - sum[j>i] Uij * xj) / Uii double sum = 0; for (int j = i + 1; j < used; ++j) sum += row.getQuick(j) * xd[j]; xd[i] = (bd[i] - sum) / row.getQuick(i); } return x; }
private DoubleMatrix1D upperTransSolve(DoubleMatrix1D b, DoubleMatrix1D x) { x.assign(b); double[] xd = ((DenseDoubleMatrix1D) x).elements(); int rows = LU.rows(); for (int i = 0; i < rows; ++i) { // Get row i SparseDoubleMatrix1D row = LU.viewRow(i); int used = (int) row.size(); // Solve for the current entry xd[i] /= row.getQuick(i); // Move this known solution over to the right hand side for the // remaining equations for (int j = i + 1; j < used; ++j) xd[j] -= row.getQuick(j) * xd[i]; } return x; }
private DoubleMatrix1D upperTransSolve(DoubleMatrix1D b, DoubleMatrix1D x) { x.assign(b); double[] xd = ((DenseDoubleMatrix1D) x).elements(); int rows = LU.rows(); for (int i = 0; i < rows; ++i) { // Get row i SparseDoubleMatrix1D row = LU.viewRow(i); int used = (int) row.size(); // Solve for the current entry xd[i] /= row.getQuick(i); // Move this known solution over to the right hand side for the // remaining equations for (int j = i + 1; j < used; ++j) xd[j] -= row.getQuick(j) * xd[i]; } return x; }
private void factor() { int n = LU.rows(); for (int i = 1; i < n; ++i) { // Get row i SparseDoubleMatrix1D rowi = LU.viewRow(i); // Drop tolerance on current row double taui = DenseDoubleAlgebra.DEFAULT.norm(rowi, Norm.Two) * tau; for (int k = 0; k < i; ++k) { // Get row k SparseDoubleMatrix1D rowk = LU.viewRow(k); if (rowk.getQuick(k) == 0) throw new RuntimeException("Zero diagonal entry on row " + (k + 1) + " during ILU process"); double LUik = rowi.getQuick(k) / rowk.getQuick(k); // Check for small elimination entry if (Math.abs(LUik) <= taui) continue; // Traverse the sparse row k, reducing row i int rowUsed = (int) rowk.size(); for (int j = k + 1; j < rowUsed; ++j) rowi.setQuick(j, rowi.getQuick(j) - LUik * rowk.getQuick(j)); // The above has overwritten LUik, so remedy that rowi.setQuick(k, LUik); } // Store back into the LU matrix, dropping as needed gather(rowi, taui, i); } // System.out.println(LU.toString()); }
private void factor() { int n = LU.rows(); for (int i = 1; i < n; ++i) { // Get row i SparseDoubleMatrix1D rowi = LU.viewRow(i); // Drop tolerance on current row double taui = DenseDoubleAlgebra.DEFAULT.norm(rowi, Norm.Two) * tau; for (int k = 0; k < i; ++k) { // Get row k SparseDoubleMatrix1D rowk = LU.viewRow(k); if (rowk.getQuick(k) == 0) throw new RuntimeException("Zero diagonal entry on row " + (k + 1) + " during ILU process"); double LUik = rowi.getQuick(k) / rowk.getQuick(k); // Check for small elimination entry if (Math.abs(LUik) <= taui) continue; // Traverse the sparse row k, reducing row i int rowUsed = (int) rowk.size(); for (int j = k + 1; j < rowUsed; ++j) rowi.setQuick(j, rowi.getQuick(j) - LUik * rowk.getQuick(j)); // The above has overwritten LUik, so remedy that rowi.setQuick(k, LUik); } // Store back into the LU matrix, dropping as needed gather(rowi, taui, i); } // System.out.println(LU.toString()); }