no.uib.cipr.matrix.sparse.ILU Maven / Gradle / Ivy

Go to download
/*
 * 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 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;
    }
}