ConjugateGradient.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.linear;

  19. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  20. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  21. import org.apache.commons.math4.legacy.exception.NullArgumentException;
  22. import org.apache.commons.math4.legacy.exception.util.ExceptionContext;

  24. /**
  25.  * <p>
  26.  * This is an implementation of the conjugate gradient method for
  27.  * {@link RealLinearOperator}. It follows closely the template by <a
  28.  * href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at
  29.  * hand is A &middot; x = b, and the residual is r = b - A &middot; x.
  30.  * </p>
  31.  * <p><b><a id="stopcrit">Default stopping criterion</a></b></p>
  32.  * <p>
  33.  * A default stopping criterion is implemented. The iterations stop when || r ||
  34.  * &le; &delta; || b ||, where b is the right-hand side vector, r the current
  35.  * estimate of the residual, and &delta; a user-specified tolerance. It should
  36.  * be noted that r is the so-called <em>updated</em> residual, which might
  37.  * differ from the true residual due to rounding-off errors (see e.g. <a
  38.  * href="#STRA2002">Strakos and Tichy, 2002</a>).
  39.  * </p>
  40.  * <p><b>Iteration count</b></p>
  41.  * <p>
  42.  * In the present context, an iteration should be understood as one evaluation
  43.  * of the matrix-vector product A &middot; x. The initialization phase therefore
  44.  * counts as one iteration.
  45.  * </p>
  46.  * <p><b><a id="context">Exception context</a></b></p>
  47.  * <p>
  48.  * Besides standard {@link DimensionMismatchException}, this class might throw
  49.  * {@link NonPositiveDefiniteOperatorException} if the linear operator or
  50.  * the preconditioner are not positive definite. In this case, the
  51.  * {@link ExceptionContext} provides some more information
  52.  * <ul>
  53.  * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
  54.  * <li>key {@code "vector"} points to the offending vector, say x, such that
  55.  * x<sup>T</sup> &middot; L &middot; x &lt; 0.</li>
  56.  * </ul>
  57.  *
  58.  * <p><b>References</b></p>
  59.  * <dl>
  60.  * <dt><a id="BARR1994">Barret et al. (1994)</a></dt>
  61.  * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra,
  62.  * V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
  63.  * <a href="http://www.netlib.org/linalg/html_templates/Templates.html"><em>
  64.  * Templates for the Solution of Linear Systems: Building Blocks for Iterative
  65.  * Methods</em></a>, SIAM</dd>
  66.  * <dt><a id="STRA2002">Strakos and Tichy (2002)</a></dt>
  67.  * <dd>Z. Strakos and P. Tichy, <a
  68.  * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf">
  69.  * <em>On error estimation in the conjugate gradient method and why it works
  70.  * in finite precision computations</em></a>, Electronic Transactions on
  71.  * Numerical Analysis 13: 56-80, 2002</dd>
  72.  * </dl>
  73.  *
  74.  * @since 3.0
  75.  */
  76. public class ConjugateGradient
  77.     extends PreconditionedIterativeLinearSolver {

  79.     /** Key for the <a href="#context">exception context</a>. */
  80.     public static final String OPERATOR = "operator";

  82.     /** Key for the <a href="#context">exception context</a>. */
  83.     public static final String VECTOR = "vector";

  85.     /**
  86.      * {@code true} if positive-definiteness of matrix and preconditioner should
  87.      * be checked.
  88.      */
  89.     private boolean check;

  91.     /** The value of &delta;, for the default stopping criterion. */
  92.     private final double delta;

  94.     /**
  95.      * Creates a new instance of this class, with <a href="#stopcrit">default
  96.      * stopping criterion</a>.
  97.      *
  98.      * @param maxIterations the maximum number of iterations
  99.      * @param delta the &delta; parameter for the default stopping criterion
  100.      * @param check {@code true} if positive definiteness of both matrix and
  101.      * preconditioner should be checked
  102.      */
  103.     public ConjugateGradient(final int maxIterations, final double delta,
  104.                              final boolean check) {
  105.         super(maxIterations);
  106.         this.delta = delta;
  107.         this.check = check;
  108.     }

  110.     /**
  111.      * Creates a new instance of this class, with <a href="#stopcrit">default
  112.      * stopping criterion</a> and custom iteration manager.
  113.      *
  114.      * @param manager the custom iteration manager
  115.      * @param delta the &delta; parameter for the default stopping criterion
  116.      * @param check {@code true} if positive definiteness of both matrix and
  117.      * preconditioner should be checked
  118.      * @throws NullArgumentException if {@code manager} is {@code null}
  119.      */
  120.     public ConjugateGradient(final IterationManager manager,
  121.                              final double delta, final boolean check)
  122.         throws NullArgumentException {
  123.         super(manager);
  124.         this.delta = delta;
  125.         this.check = check;
  126.     }

  128.     /**
  129.      * Returns {@code true} if positive-definiteness should be checked for both
  130.      * matrix and preconditioner.
  131.      *
  132.      * @return {@code true} if the tests are to be performed
  133.      */
  134.     public final boolean getCheck() {
  135.         return check;
  136.     }

  138.     /**
  139.      * {@inheritDoc}
  140.      *
  141.      * @throws NonPositiveDefiniteOperatorException if {@code a} or {@code m} is
  142.      * not positive definite
  143.      */
  144.     @Override
  145.     public RealVector solveInPlace(final RealLinearOperator a,
  146.                                    final RealLinearOperator m,
  147.                                    final RealVector b,
  148.                                    final RealVector x0)
  149.         throws NullArgumentException, NonPositiveDefiniteOperatorException,
  150.         NonSquareOperatorException, DimensionMismatchException,
  151.         MaxCountExceededException {
  152.         checkParameters(a, m, b, x0);
  153.         final IterationManager manager = getIterationManager();
  154.         // Initialization of default stopping criterion
  155.         manager.resetIterationCount();
  156.         final double rmax = delta * b.getNorm();
  157.         final RealVector bro = RealVector.unmodifiableRealVector(b);

  159.         // Initialization phase counts as one iteration.
  160.         manager.incrementIterationCount();
  161.         // p and x are constructed as copies of x0, since presumably, the type
  162.         // of x is optimized for the calculation of the matrix-vector product
  163.         // A.x.
  164.         final RealVector x = x0;
  165.         final RealVector xro = RealVector.unmodifiableRealVector(x);
  166.         final RealVector p = x.copy();
  167.         RealVector q = a.operate(p);

  169.         final RealVector r = b.combine(1, -1, q);
  170.         final RealVector rro = RealVector.unmodifiableRealVector(r);
  171.         double rnorm = r.getNorm();
  172.         RealVector z;
  173.         if (m == null) {
  174.             z = r;
  175.         } else {
  176.             z = null;
  177.         }
  178.         IterativeLinearSolverEvent evt;
  179.         evt = new DefaultIterativeLinearSolverEvent(this,
  180.             manager.getIterations(), xro, bro, rro, rnorm);
  181.         manager.fireInitializationEvent(evt);
  182.         if (rnorm <= rmax) {
  183.             manager.fireTerminationEvent(evt);
  184.             return x;
  185.         }
  186.         double rhoPrev = 0.;
  187.         while (true) {
  188.             manager.incrementIterationCount();
  189.             evt = new DefaultIterativeLinearSolverEvent(this,
  190.                 manager.getIterations(), xro, bro, rro, rnorm);
  191.             manager.fireIterationStartedEvent(evt);
  192.             if (m != null) {
  193.                 z = m.operate(r);
  194.             }
  195.             final double rhoNext = r.dotProduct(z);
  196.             if (check && rhoNext <= 0) {
  197.                 final NonPositiveDefiniteOperatorException e;
  198.                 e = new NonPositiveDefiniteOperatorException();
  199.                 final ExceptionContext context = e.getContext();
  200.                 context.setValue(OPERATOR, m);
  201.                 context.setValue(VECTOR, r);
  202.                 throw e;
  203.             }
  204.             if (manager.getIterations() == 2) {
  205.                 p.setSubVector(0, z);
  206.             } else {
  207.                 p.combineToSelf(rhoNext / rhoPrev, 1., z);
  208.             }
  209.             q = a.operate(p);
  210.             final double pq = p.dotProduct(q);
  211.             if (check && pq <= 0) {
  212.                 final NonPositiveDefiniteOperatorException e;
  213.                 e = new NonPositiveDefiniteOperatorException();
  214.                 final ExceptionContext context = e.getContext();
  215.                 context.setValue(OPERATOR, a);
  216.                 context.setValue(VECTOR, p);
  217.                 throw e;
  218.             }
  219.             final double alpha = rhoNext / pq;
  220.             x.combineToSelf(1., alpha, p);
  221.             r.combineToSelf(1., -alpha, q);
  222.             rhoPrev = rhoNext;
  223.             rnorm = r.getNorm();
  224.             evt = new DefaultIterativeLinearSolverEvent(this,
  225.                 manager.getIterations(), xro, bro, rro, rnorm);
  226.             manager.fireIterationPerformedEvent(evt);
  227.             if (rnorm <= rmax) {
  228.                 manager.fireTerminationEvent(evt);
  229.                 return x;
  230.             }
  231.         }
  232.     }
  233. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT