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 · x = b, and the residual is r = b - A · 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 * ≤ δ || b ||, where b is the right-hand side vector, r the current
036 * estimate of the residual, and δ 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 · 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> · L · 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 δ, 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 δ 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 δ 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 }