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

Go to download
Show more of this group Show more artifacts with this name
Show all versions of mt-java Show documentation
Matrix data structures, linear solvers, least squares methods, eigenvalue, and singular value decompositions. For larger random dense matrices (above ~ 350 x 350) matrix-matrix multiplication C = A.B is about 50% faster than MTJ.
The newest version!
/*
 * 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
 */
/*
 * Derived from public domain software at http://www.netlib.org/templates
 */

package no.uib.cipr.matrix.sparse;

import no.uib.cipr.matrix.DenseMatrix;
import no.uib.cipr.matrix.DenseVector;
import no.uib.cipr.matrix.GivensRotation;
import no.uib.cipr.matrix.Matrix;
import no.uib.cipr.matrix.UpperTriangDenseMatrix;
import no.uib.cipr.matrix.Vector;
import no.uib.cipr.matrix.Vector.Norm;

/**
 * GMRES solver. GMRES solves the unsymmetric linear system Ax = b
 * using the Generalized Minimum Residual method. The GMRES iteration is
 * restarted after a given number of iterations. By default it is restarted
 * after 30 iterations.
 * 
 * @author Templates
 */
public class GMRES extends AbstractIterativeSolver {

    /**
     * After this many iterations, the GMRES will be restarted.
     */
    private int restart;

    /**
     * Vectors for use in the iterative solution process
     */
    private Vector w, u, r;

    /**
     * Vectors spanning the subspace
     */
    private Vector[] v;

    /**
     * Restart vector
     */
    private DenseVector s;

    /**
     * Hessenberg matrix
     */
    private DenseMatrix H;

    /**
     * Givens rotations for the QR factorization
     */
    private GivensRotation[] rotation;

    /**
     * Constructor for GMRES. Uses the given vector as template for creating
     * scratch vectors. Typically, the solution or the right hand side vector
     * can be passed, and the template is not modified. The iteration is
     * restarted every 30 iterations.
     * 
     * @param template
     *            Vector to use as template for the work vectors needed in the
     *            solution process
     */
    public GMRES(Vector template) {
        this(template, 30);
    }

    /**
     * Constructor for GMRES. Uses the given vector as template for creating
     * scratch vectors. Typically, the solution or the right hand side vector
     * can be passed, and the template is not modified.
     * 
     * @param template
     *            Vector to use as template for the work vectors needed in the
     *            solution process
     * @param restart
     *            GMRES iteration is restarted after this number of iterations
     */
    public GMRES(Vector template, int restart) {
        w = template.copy();
        u = template.copy();
        r = template.copy();
        setRestart(restart);
    }

    /**
     * Sets the restart parameter.
     * 
     * @param restart
     *            GMRES iteration is restarted after this number of iterations
     */
    public void setRestart(int restart) {
        this.restart = restart;
        if (restart <= 0)
            throw new IllegalArgumentException(
                    "restart must be a positive integer");

        s = new DenseVector(restart + 1);
        H = new DenseMatrix(restart + 1, restart);
        rotation = new GivensRotation[restart + 1];

        v = new Vector[restart + 1];
        for (int i = 0; i < v.length; ++i)
            v[i] = r.copy().zero();
    }

    public Vector solve(Matrix A, Vector b, Vector x)
            throws IterativeSolverNotConvergedException {
        checkSizes(A, b, x);

        A.multAdd(-1, x, u.set(b));
        M.apply(u, r);
        double normr = r.norm(Norm.Two);
        M.apply(b, u);

        // Outer iteration
        for (iter.setFirst(); !iter.converged(r, x); iter.next()) {

            v[0].set(1 / normr, r);
            s.zero().set(0, normr);
            int i = 0;

            // Inner iteration
            for (; i < restart && !iter.converged(Math.abs(s.get(i))); i++, iter
                    .next()) {
                A.mult(v[i], u);
                M.apply(u, w);

                for (int k = 0; k <= i; k++) {
                    H.set(k, i, w.dot(v[k]));
                    w.add(-H.get(k, i), v[k]);
                }
                H.set(i + 1, i, w.norm(Norm.Two));
                v[i + 1].set(1. / H.get(i + 1, i), w);

                // QR factorization of H using Givens rotations
                for (int k = 0; k < i; ++k)
                    rotation[k].apply(H, i, k, k + 1);

                rotation[i] = new GivensRotation(H.get(i, i), H.get(i + 1, i));
                rotation[i].apply(H, i, i, i + 1);
                rotation[i].apply(s, i, i + 1);
            }

            // Update solution in current subspace
            new UpperTriangDenseMatrix(H, i, false).solve(s, s);
            for (int j = 0; j < i; j++)
                x.add(s.get(j), v[j]);

            A.multAdd(-1, x, u.set(b));
            M.apply(u, r);
            normr = r.norm(Norm.Two);
        }

        return x;
    }
}