/* * Copyright (C) 2003-2006 Bjørn-Ove Heimsund * * This file is part of MTJ. * * This library is free software; you can redistribute it and/or modify it * under the terms of the GNU Lesser General Public License as published by the * Free Software Foundation; either version 2.1 of the License, or (at your * option) any later version. * * This library is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this library; if not, write to the Free Software Foundation, * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ package no.uib.cipr.matrix.sparse; import no.uib.cipr.matrix.DenseVector; import no.uib.cipr.matrix.Matrix; import no.uib.cipr.matrix.Vector; /** * ILU(0) preconditioner using a compressed row matrix as internal storage */ public class ILU implements Preconditioner { /** * Factorisation matrix */ private final CompRowMatrix LU; /** * The L and U factors */ private Matrix L, U; /** * Temporary vector for solving the factorised system */ private final Vector y; /** * Sets up the ILU preconditioner * * @param LU * Matrix to use internally. For best performance, its non-zero * pattern must conform to that of the system matrix */ public ILU(CompRowMatrix LU) { if (!LU.isSquare()) throw new IllegalArgumentException( "ILU only applies to square matrices"); this.LU = LU; int n = LU.numRows(); y = new DenseVector(n); } public Vector apply(Vector b, Vector x) { // Ly = b, y = L\b L.solve(b, y); // Ux = L\b = y return U.solve(y, x); } public Vector transApply(Vector b, Vector x) { // U'y = b, y = U'\b U.transSolve(b, y); // L'x = U'\b = y return L.transSolve(y, x); } public void setMatrix(Matrix A) { LU.set(A); factor(); } private void factor() { int n = LU.numRows(); // Internal CRS matrix storage int[] colind = LU.getColumnIndices(); int[] rowptr = LU.getRowPointers(); double[] data = LU.getData(); // Find the indices to the diagonal entries int[] diagind = findDiagonalIndices(n, colind, rowptr); // Go down along the main diagonal for (int k = 1; k < n; ++k) for (int i = rowptr[k]; i < diagind[k]; ++i) { // Get the current diagonal entry int index = colind[i]; double LUii = data[diagind[index]]; if (LUii == 0) throw new RuntimeException("Zero pivot encountered on row " + (i + 1) + " during ILU process"); // Elimination factor double LUki = (data[i] /= LUii); // Traverse the sparse row i, reducing on row k for (int j = diagind[index] + 1, l = rowptr[k] + 1; j < rowptr[index + 1]; ++j) { while (l < rowptr[k + 1] && colind[l] < colind[j]) l++; if (colind[l] == colind[j]) data[l] -= LUki * data[j]; } } L = new UnitLowerCompRowMatrix(LU, diagind); U = new UpperCompRowMatrix(LU, diagind); } private static int[] findDiagonalIndices(int m, int[] colind, int[] rowptr) { int[] diagind = new int[m]; for (int k = 0; k < m; ++k) { diagind[k] = Arrays.binarySearch(colind, k, rowptr[k], rowptr[k + 1]); if (diagind[k] < 0) throw new RuntimeException("Missing diagonal entry on row " + (k + 1)); } return diagind; } }