amg.java

另一个功能更强大的矩阵运算软件开源代码

/*
 * 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 java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.Set;

import no.uib.cipr.matrix.DenseLU;
import no.uib.cipr.matrix.DenseMatrix;
import no.uib.cipr.matrix.DenseVector;
import no.uib.cipr.matrix.Matrix;
import no.uib.cipr.matrix.Vector;

/**
 * Algebraic multigrid preconditioner. Uses the smoothed aggregation method
 * described by Vanek, Mandel, and Brezina (1996).
 */
public class AMG implements Preconditioner {

    /**
     * Relaxations at each level
     */
    private SSOR[] preM, postM;

    /**
     * The number of levels
     */
    private int m;

    /**
     * System matrix at each level, except at the coarsest
     */
    private CompRowMatrix[] A;

    /**
     * LU factorization at the coarsest level
     */
    private DenseLU lu;

    /**
     * Solution, right-hand side, and residual vectors at each level
     */
    private DenseVector[] u, f, r;

    /**
     * Interpolation operators going to a finer mesh
     */
    private CompColMatrix[] I;

    /**
     * Smallest matrix size before terminating the AMG setup phase. Matrices
     * smaller than this will be solved by a direct solver
     */
    private final int min;

    /**
     * Number of times to perform the pre- and post-smoothings
     */
    private final int nu1, nu2;

    /**
     * Determines cycle type. gamma=1 is V, gamma=2 is W
     */
    private final int gamma;

    /**
     * Overrelaxation parameters in the pre- and post-smoothings, and with the
     * possibility of distinct values in the forward and reverse sweeps
     */
    private final double omegaPreF, omegaPreR, omegaPostF, omegaPostR;

    /**
     * Perform a reverse (backwards) smoothing sweep
     */
    private final boolean reverse;

    /**
     * Jacobi damping parameter, between zero and one. If it equals zero, the
     * method reduces to the standard aggregate multigrid method
     */
    private final double omega;

    /**
     * Operating in transpose mode?
     */
    private boolean transpose;

    /**
     * Sets up the algebraic multigrid preconditioner
     * 
     * @param omegaPreF
     *            Overrelaxation parameter in the forward sweep of the
     *            pre-smoothing
     * @param omegaPreR
     *            Overrelaxation parameter in the backwards sweep of the
     *            pre-smoothing
     * @param omegaPostF
     *            Overrelaxation parameter in the forward sweep of the
     *            post-smoothing
     * @param omegaPostR
     *            Overrelaxation parameter in the backwards sweep of the
     *            post-smoothing
     * @param nu1
     *            Number of pre-relaxations to perform
     * @param nu2
     *            Number of post-relaxations to perform
     * @param gamma
     *            Number of times to go to a coarser level
     * @param min
     *            Smallest matrix size before using a direct solver
     * @param omega
     *            Jacobi damping parameter, between zero and one. If it equals
     *            zero, the method reduces to the standard aggregate multigrid
     *            method
     */
    public AMG(double omegaPreF, double omegaPreR, double omegaPostF,
            double omegaPostR, int nu1, int nu2, int gamma, int min,
            double omega) {
        this.omegaPreF = omegaPreF;
        this.omegaPreR = omegaPreR;
        this.omegaPostF = omegaPostF;
        this.omegaPostR = omegaPostR;

        reverse = true;

        this.nu1 = nu1;
        this.nu2 = nu2;
        this.gamma = gamma;
        this.min = min;

        this.omega = omega;
    }

    /**
     * Sets up the algebraic multigrid preconditioner. Uses an SOR method,
     * without the backward sweep in SSOR
     * 
     * @param omegaPre
     *            Overrelaxation parameter in the pre-smoothing
     * @param omegaPost
     *            Overrelaxation parameter in the post-smoothing
     * @param nu1
     *            Number of pre-relaxations to perform
     * @param nu2
     *            Number of post-relaxations to perform
     * @param gamma
     *            Number of times to go to a coarser level
     * @param min
     *            Smallest matrix size before using a direct solver
     * @param omega
     *            Jacobi damping parameter, between zero and one. If it equals
     *            zero, the method reduces to the standard aggregate multigrid
     *            method
     */
    public AMG(double omegaPre, double omegaPost, int nu1, int nu2, int gamma,
            int min, double omega) {
        this.omegaPreF = omegaPre;
        this.omegaPreR = omegaPre;
        this.omegaPostF = omegaPost;
        this.omegaPostR = omegaPost;

        reverse = false;

        this.nu1 = nu1;
        this.nu2 = nu2;
        this.gamma = gamma;
        this.min = min;

        this.omega = omega;
    }

    /**
     * Sets up the algebraic multigrid preconditioner using some default
     * parameters. In the presmoothing, <code>omegaF=1</code> and
     * <code>omegaR=1.85</code>, while in the postsmoothing,
     * <code>omegaF=1.85</code> and <code>omegaR=1</code>. Sets
     * <code>nu1=nu2=gamma=1</code>, has a smallest matrix size of 40, and
     * sets <code>omega=2/3</code>.
     */
    public AMG() {
        this(1, 1.85, 1.85, 1, 1, 1, 1, 40, 2. / 3);
    }

    public Vector apply(Vector b, Vector x) {
        u[0].set(x);
        f[0].set(b);

        transpose = false;
        cycle(0);

        return x.set(u[0]);
    }

    public Vector transApply(Vector b, Vector x) {
        u[0].set(x);
        f[0].set(b);

        transpose = true;
        cycle(0);

        return x.set(u[0]);
    }

    public void setMatrix(Matrix A) {
        List<CompRowMatrix> Al = new LinkedList<CompRowMatrix>();
        List<CompColMatrix> Il = new LinkedList<CompColMatrix>();

        Al.add(new CompRowMatrix(A));

        for (int k = 0; Al.get(k).numRows() > min; ++k) {

            CompRowMatrix Af = Al.get(k);

            double eps = 0.08 * Math.pow(0.5, k);

            // Create the aggregates
            Aggregator aggregator = new Aggregator(Af, eps);

            // If no aggregates were created, no interpolation operator will be
            // created, and the setup phase stops
            if (aggregator.getAggregates().size() == 0)
                break;

            // Create an interpolation operator using smoothing. This also
            // creates the Galerkin operator
            Interpolator sa = new Interpolator(aggregator, Af, omega);

            Al.add(sa.getGalerkinOperator());
            Il.add(sa.getInterpolationOperator());
        }

        // Copy to array storage
        m = Al.size();
        if (m == 0)
            throw new RuntimeException("Matrix too small for AMG");

        I = new CompColMatrix[m - 1];
        this.A = new CompRowMatrix[m - 1];

        Il.toArray(I);
        for (int i = 0; i < Al.size() - 1; ++i)
            this.A[i] = Al.get(i);

        // Create a LU decomposition of the smallest Galerkin matrix
        DenseMatrix Ac = new DenseMatrix(Al.get(Al.size() - 1));
        lu = new DenseLU(Ac.numRows(), Ac.numColumns());
        lu.factor(Ac);

        // Allocate vectors at each level
        u = new DenseVector[m];
        f = new DenseVector[m];
        r = new DenseVector[m];
        for (int k = 0; k < m; ++k) {
            int n = Al.get(k).numRows();
            u[k] = new DenseVector(n);
            f[k] = new DenseVector(n);
            r[k] = new DenseVector(n);
        }

        // Set up the SSOR relaxation schemes
        preM = new SSOR[m - 1];
        postM = new SSOR[m - 1];
        for (int k = 0; k < m - 1; ++k) {
            CompRowMatrix Ak = this.A[k];
            preM[k] = new SSOR(Ak, reverse, omegaPreF, omegaPreR);
            postM[k] = new SSOR(Ak, reverse, omegaPostF, omegaPostR);
            preM[k].setMatrix(Ak);
            postM[k].setMatrix(Ak);
        }
    }

    /**
     * Performs a multigrid cycle
     * 
     * @param k
     *            Level to cycle at. Start by calling <code>cycle(0)</code>
     */
    private void cycle(int k) {
        if (k == m - 1)
            directSolve();
        else {

            // Presmoothings
            preRelax(k);

            u[k + 1].zero();

            // Compute the residual
            A[k].multAdd(-1, u[k], r[k].set(f[k]));

            // Restrict to the next coarser level
            I[k].transMult(r[k], f[k + 1]);

            // Recurse to next level
            for (int i = 0; i < gamma; ++i)
                cycle(k + 1);

            // Add residual correction by prolongation
            I[k].multAdd(u[k + 1], u[k]);

            // Postsmoothings
            postRelax(k);
        }
    }

    /**
     * Solves directly at the coarsest level
     */
    private void directSolve() {
        int k = m - 1;
        u[k].set(f[k]);
        DenseMatrix U = new DenseMatrix(u[k], false);

        if (transpose)
            lu.transSolve(U);
        else
            lu.solve(U);
    }

    /**
     * Applies the relaxation scheme at the given level
     * 
     * @param k
     *            Multigrid level
     */
    private void preRelax(int k) {
        for (int i = 0; i < nu1; ++i)
            if (transpose)
                preM[k].transApply(f[k], u[k]);
            else
                preM[k].apply(f[k], u[k]);
    }

    /**
     * Applies the relaxation scheme at the given level
     * 
     * @param k
     *            Multigrid level
     */
    private void postRelax(int k) {
        for (int i = 0; i < nu2; ++i)
            if (transpose)
                postM[k].transApply(f[k], u[k]);
            else
                postM[k].apply(f[k], u[k]);
    }

    /**
     * Creates aggregates. These are disjoint sets, each of which represents one
     * node at a coarser mesh by aggregating together a set of fine nodes
     */
    private static class Aggregator {

        /**
         * The aggregates
         */
        private List<Set<Integer>> C;

        /**
         * Diagonal indices into the sparse matrix
         */
        private int[] diagind;

        /**
         * The strongly coupled node neighborhood of a given node
         */
        private List<Set<Integer>> N;

        /**
         * Creates the aggregates
         * 
         * @param A
         *            Sparse matrix
         * @param eps
         *            Tolerance for selecting the strongly coupled node
         *            neighborhoods. Between zero and one.
         */
        public Aggregator(CompRowMatrix A, double eps) {

            diagind = findDiagonalIndices(A);
            N = findNodeNeighborhood(A, diagind, eps);

            /*
             * Initialization. Remove isolated nodes from the aggregates
             */

            boolean[] R = createInitialR(A);

            /*
             * Startup aggregation. Use disjoint strongly coupled neighborhoods
             * as the initial aggregate approximation
             */

            C = createInitialAggregates(N, R);

            /*
             * Enlargment of the aggregates. Add nodes to each aggregate based
             * on how strongly connected the nodes are to a given aggregate
             */

            C = enlargeAggregates(C, N, R);

            /*
             * Handling of the remenants. Put all remaining unallocated nodes
             * into new aggregates defined by the intersection of N and R
             */

            C = createFinalAggregates(C, N, R);
        }

        /**
         * Gets the aggregates
         */
        public List<Set<Integer>> getAggregates() {
            return C;
        }

        /**
         * Returns the matrix diagonal indices. This is a by-product of the
         * aggregation
         */
        public int[] getDiagonalIndices() {
            return diagind;
        }

        /**
         * Returns the strongly coupled node neighborhoods of a given node. This
         * is a by-product of the aggregation
         */
        public List<Set<Integer>> getNodeNeighborhoods() {
            return N;
        }

        /**
         * Finds the diagonal indices of the matrix
         */
        private int[] findDiagonalIndices(CompRowMatrix A) {
            int[] rowptr = A.getRowPointers();
            int[] colind = A.getColumnIndices();

            int[] diagind = new int[A.numRows()];