    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.linear;
    018
    019    import org.apache.commons.math3.exception.DimensionMismatchException;
    020    import org.apache.commons.math3.exception.MaxCountExceededException;
    021    import org.apache.commons.math3.exception.NullArgumentException;
    022    import org.apache.commons.math3.exception.util.ExceptionContext;
    023    import org.apache.commons.math3.util.IterationManager;
    024
    025    /**
    026     * <p>
    027     * This is an implementation of the conjugate gradient method for
    028     * {@link RealLinearOperator}. It follows closely the template by <a
    029     * href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at
    030     * hand is A &middot; x = b, and the residual is r = b - A &middot; x.
    031     * </p>
    032     * <h3><a id="stopcrit">Default stopping criterion</a></h3>
    033     * <p>
    034     * A default stopping criterion is implemented. The iterations stop when || r ||
    035     * &le; &delta; || b ||, where b is the right-hand side vector, r the current
    036     * estimate of the residual, and &delta; a user-specified tolerance. It should
    037     * be noted that r is the so-called <em>updated</em> residual, which might
    038     * differ from the true residual due to rounding-off errors (see e.g. <a
    039     * href="#STRA2002">Strakos and Tichy, 2002</a>).
    040     * </p>
    041     * <h3>Iteration count</h3>
    042     * <p>
    043     * In the present context, an iteration should be understood as one evaluation
    044     * of the matrix-vector product A &middot; x. The initialization phase therefore
    045     * counts as one iteration.
    046     * </p>
    047     * <h3><a id="context">Exception context</a></h3>
    048     * <p>
    049     * Besides standard {@link DimensionMismatchException}, this class might throw
    050     * {@link NonPositiveDefiniteOperatorException} if the linear operator or
    051     * the preconditioner are not positive definite. In this case, the
    052     * {@link ExceptionContext} provides some more information
    053     * <ul>
    054     * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
    055     * <li>key {@code "vector"} points to the offending vector, say x, such that
    056     * x<sup>T</sup> &middot; L &middot; x < 0.</li>
    057     * </ul>
    058     * </p>
    059     * <h3>References</h3>
    060     * <dl>
    061     * <dt><a id="BARR1994">Barret et al. (1994)</a></dt>
    062     * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra,
    063     * V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
    064     * <a href="http://www.netlib.org/linalg/html_templates/Templates.html"><em>
    065     * Templates for the Solution of Linear Systems: Building Blocks for Iterative
    066     * Methods</em></a>, SIAM</dd>
    067     * <dt><a id="STRA2002">Strakos and Tichy (2002)
    068     * <dt>
    069     * <dd>Z. Strakos and P. Tichy, <a
    070     * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf">
    071     * <em>On error estimation in the conjugate gradient method and why it works
    072     * in finite precision computations</em></a>, Electronic Transactions on
    073     * Numerical Analysis 13: 56-80, 2002</dd>
    074     * </dl>
    075     *
    076     * @version $Id: ConjugateGradient.java 1416643 2012-12-03 19:37:14Z tn $
    077     * @since 3.0
    078     */
    079    public class ConjugateGradient
    080        extends PreconditionedIterativeLinearSolver {
    081
    082        /** Key for the <a href="#context">exception context</a>. */
    083        public static final String OPERATOR = "operator";
    084
    085        /** Key for the <a href="#context">exception context</a>. */
    086        public static final String VECTOR = "vector";
    087
    088        /**
    089         * {@code true} if positive-definiteness of matrix and preconditioner should
    090         * be checked.
    091         */
    092        private boolean check;
    093
    094        /** The value of &delta;, for the default stopping criterion. */
    095        private final double delta;
    096
    097        /**
    098         * Creates a new instance of this class, with <a href="#stopcrit">default
    099         * stopping criterion</a>.
    100         *
    101         * @param maxIterations the maximum number of iterations
    102         * @param delta the &delta; parameter for the default stopping criterion
    103         * @param check {@code true} if positive definiteness of both matrix and
    104         * preconditioner should be checked
    105         */
    106        public ConjugateGradient(final int maxIterations, final double delta,
    107                                 final boolean check) {
    108            super(maxIterations);
    109            this.delta = delta;
    110            this.check = check;
    111        }
    112
    113        /**
    114         * Creates a new instance of this class, with <a href="#stopcrit">default
    115         * stopping criterion</a> and custom iteration manager.
    116         *
    117         * @param manager the custom iteration manager
    118         * @param delta the &delta; parameter for the default stopping criterion
    119         * @param check {@code true} if positive definiteness of both matrix and
    120         * preconditioner should be checked
    121         * @throws NullArgumentException if {@code manager} is {@code null}
    122         */
    123        public ConjugateGradient(final IterationManager manager,
    124                                 final double delta, final boolean check)
    125            throws NullArgumentException {
    126            super(manager);
    127            this.delta = delta;
    128            this.check = check;
    129        }
    130
    131        /**
    132         * Returns {@code true} if positive-definiteness should be checked for both
    133         * matrix and preconditioner.
    134         *
    135         * @return {@code true} if the tests are to be performed
    136         */
    137        public final boolean getCheck() {
    138            return check;
    139        }
    140
    141        /**
    142         * {@inheritDoc}
    143         *
    144         * @throws NonPositiveDefiniteOperatorException if {@code a} or {@code m} is
    145         * not positive definite
    146         */
    147        @Override
    148        public RealVector solveInPlace(final RealLinearOperator a,
    149                                       final RealLinearOperator m,
    150                                       final RealVector b,
    151                                       final RealVector x0)
    152            throws NullArgumentException, NonPositiveDefiniteOperatorException,
    153            NonSquareOperatorException, DimensionMismatchException,
    154            MaxCountExceededException, NonPositiveDefiniteOperatorException {
    155            checkParameters(a, m, b, x0);
    156            final IterationManager manager = getIterationManager();
    157            // Initialization of default stopping criterion
    158            manager.resetIterationCount();
    159            final double rmax = delta * b.getNorm();
    160            final RealVector bro = RealVector.unmodifiableRealVector(b);
    161
    162            // Initialization phase counts as one iteration.
    163            manager.incrementIterationCount();
    164            // p and x are constructed as copies of x0, since presumably, the type
    165            // of x is optimized for the calculation of the matrix-vector product
    166            // A.x.
    167            final RealVector x = x0;
    168            final RealVector xro = RealVector.unmodifiableRealVector(x);
    169            final RealVector p = x.copy();
    170            RealVector q = a.operate(p);
    171
    172            final RealVector r = b.combine(1, -1, q);
    173            final RealVector rro = RealVector.unmodifiableRealVector(r);
    174            double rnorm = r.getNorm();
    175            RealVector z;
    176            if (m == null) {
    177                z = r;
    178            } else {
    179                z = null;
    180            }
    181            IterativeLinearSolverEvent evt;
    182            evt = new DefaultIterativeLinearSolverEvent(this,
    183                manager.getIterations(), xro, bro, rro, rnorm);
    184            manager.fireInitializationEvent(evt);
    185            if (rnorm <= rmax) {
    186                manager.fireTerminationEvent(evt);
    187                return x;
    188            }
    189            double rhoPrev = 0.;
    190            while (true) {
    191                manager.incrementIterationCount();
    192                evt = new DefaultIterativeLinearSolverEvent(this,
    193                    manager.getIterations(), xro, bro, rro, rnorm);
    194                manager.fireIterationStartedEvent(evt);
    195                if (m != null) {
    196                    z = m.operate(r);
    197                }
    198                final double rhoNext = r.dotProduct(z);
    199                if (check && (rhoNext <= 0.)) {
    200                    final NonPositiveDefiniteOperatorException e;
    201                    e = new NonPositiveDefiniteOperatorException();
    202                    final ExceptionContext context = e.getContext();
    203                    context.setValue(OPERATOR, m);
    204                    context.setValue(VECTOR, r);
    205                    throw e;
    206                }
    207                if (manager.getIterations() == 2) {
    208                    p.setSubVector(0, z);
    209                } else {
    210                    p.combineToSelf(rhoNext / rhoPrev, 1., z);
    211                }
    212                q = a.operate(p);
    213                final double pq = p.dotProduct(q);
    214                if (check && (pq <= 0.)) {
    215                    final NonPositiveDefiniteOperatorException e;
    216                    e = new NonPositiveDefiniteOperatorException();
    217                    final ExceptionContext context = e.getContext();
    218                    context.setValue(OPERATOR, a);
    219                    context.setValue(VECTOR, p);
    220                    throw e;
    221                }
    222                final double alpha = rhoNext / pq;
    223                x.combineToSelf(1., alpha, p);
    224                r.combineToSelf(1., -alpha, q);
    225                rhoPrev = rhoNext;
    226                rnorm = r.getNorm();
    227                evt = new DefaultIterativeLinearSolverEvent(this,
    228                    manager.getIterations(), xro, bro, rro, rnorm);
    229                manager.fireIterationPerformedEvent(evt);
    230                if (rnorm <= rmax) {
    231                    manager.fireTerminationEvent(evt);
    232                    return x;
    233                }
    234            }
    235        }
    236    }