/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.linear;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.util.ExceptionContext;
import org.apache.commons.math3.util.IterationManager;
/**
* <p>
* This is an implementation of the conjugate gradient method for
* {@link RealLinearOperator}. It follows closely the template by <a
* href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at
* hand is A · x = b, and the residual is r = b - A · x.
* </p>
* <h3><a id="stopcrit">Default stopping criterion</a></h3>
* <p>
* A default stopping criterion is implemented. The iterations stop when || r ||
* ≤ δ || b ||, where b is the right-hand side vector, r the current
* estimate of the residual, and δ a user-specified tolerance. It should
* be noted that r is the so-called <em>updated</em> residual, which might
* differ from the true residual due to rounding-off errors (see e.g. <a
* href="#STRA2002">Strakos and Tichy, 2002</a>).
* </p>
* <h3>Iteration count</h3>
* <p>
* In the present context, an iteration should be understood as one evaluation
* of the matrix-vector product A · x. The initialization phase therefore
* counts as one iteration.
* </p>
* <h3><a id="context">Exception context</a></h3>
* <p>
* Besides standard {@link DimensionMismatchException}, this class might throw
* {@link NonPositiveDefiniteOperatorException} if the linear operator or
* the preconditioner are not positive definite. In this case, the
* {@link ExceptionContext} provides some more information
* <ul>
* <li>key {@code "operator"} points to the offending linear operator, say L,</li>
* <li>key {@code "vector"} points to the offending vector, say x, such that
* x<sup>T</sup> · L · x < 0.</li>
* </ul>
* </p>
* <h3>References</h3>
* <dl>
* <dt><a id="BARR1994">Barret et al. (1994)</a></dt>
* <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra,
* V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
* <a href="http://www.netlib.org/linalg/html_templates/Templates.html"><em>
* Templates for the Solution of Linear Systems: Building Blocks for Iterative
* Methods</em></a>, SIAM</dd>
* <dt><a id="STRA2002">Strakos and Tichy (2002)
* <dt>
* <dd>Z. Strakos and P. Tichy, <a
* href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf">
* <em>On error estimation in the conjugate gradient method and why it works
* in finite precision computations</em></a>, Electronic Transactions on
* Numerical Analysis 13: 56-80, 2002</dd>
* </dl>
*
* @since 3.0
*/
public class ConjugateGradient
extends PreconditionedIterativeLinearSolver {
/** Key for the <a href="#context">exception context</a>. */
public static final String OPERATOR = "operator";
/** Key for the <a href="#context">exception context</a>. */
public static final String VECTOR = "vector";
/**
* {@code true} if positive-definiteness of matrix and preconditioner should
* be checked.
*/
private boolean check;
/** The value of δ, for the default stopping criterion. */
private final double delta;
/**
* Creates a new instance of this class, with <a href="#stopcrit">default
* stopping criterion</a>.
*
* @param maxIterations the maximum number of iterations
* @param delta the δ parameter for the default stopping criterion
* @param check {@code true} if positive definiteness of both matrix and
* preconditioner should be checked
*/
public ConjugateGradient(final int maxIterations, final double delta,
final boolean check) {
super(maxIterations);
this.delta = delta;
this.check = check;
}
/**
* Creates a new instance of this class, with <a href="#stopcrit">default
* stopping criterion</a> and custom iteration manager.
*
* @param manager the custom iteration manager
* @param delta the δ parameter for the default stopping criterion
* @param check {@code true} if positive definiteness of both matrix and
* preconditioner should be checked
* @throws NullArgumentException if {@code manager} is {@code null}
*/
public ConjugateGradient(final IterationManager manager,
final double delta, final boolean check)
throws NullArgumentException {
super(manager);
this.delta = delta;
this.check = check;
}
/**
* Returns {@code true} if positive-definiteness should be checked for both
* matrix and preconditioner.
*
* @return {@code true} if the tests are to be performed
*/
public final boolean getCheck() {
return check;
}
/**
* {@inheritDoc}
*
* @throws NonPositiveDefiniteOperatorException if {@code a} or {@code m} is
* not positive definite
*/
@Override
public RealVector solveInPlace(final RealLinearOperator a,
final RealLinearOperator m,
final RealVector b,
final RealVector x0)
throws NullArgumentException, NonPositiveDefiniteOperatorException,
NonSquareOperatorException, DimensionMismatchException,
MaxCountExceededException {
checkParameters(a, m, b, x0);
final IterationManager manager = getIterationManager();
// Initialization of default stopping criterion
manager.resetIterationCount();
final double rmax = delta * b.getNorm();
final RealVector bro = RealVector.unmodifiableRealVector(b);
// Initialization phase counts as one iteration.
manager.incrementIterationCount();
// p and x are constructed as copies of x0, since presumably, the type
// of x is optimized for the calculation of the matrix-vector product
// A.x.
final RealVector x = x0;
final RealVector xro = RealVector.unmodifiableRealVector(x);
final RealVector p = x.copy();
RealVector q = a.operate(p);
final RealVector r = b.combine(1, -1, q);
final RealVector rro = RealVector.unmodifiableRealVector(r);
double rnorm = r.getNorm();
RealVector z;
if (m == null) {
z = r;
} else {
z = null;
}
IterativeLinearSolverEvent evt;
evt = new DefaultIterativeLinearSolverEvent(this,
manager.getIterations(), xro, bro, rro, rnorm);
manager.fireInitializationEvent(evt);
if (rnorm <= rmax) {
manager.fireTerminationEvent(evt);
return x;
}
double rhoPrev = 0.;
while (true) {
manager.incrementIterationCount();
evt = new DefaultIterativeLinearSolverEvent(this,
manager.getIterations(), xro, bro, rro, rnorm);
manager.fireIterationStartedEvent(evt);
if (m != null) {
z = m.operate(r);
}
final double rhoNext = r.dotProduct(z);
if (check && (rhoNext <= 0.)) {
final NonPositiveDefiniteOperatorException e;
e = new NonPositiveDefiniteOperatorException();
final ExceptionContext context = e.getContext();
context.setValue(OPERATOR, m);
context.setValue(VECTOR, r);
throw e;
}
if (manager.getIterations() == 2) {
p.setSubVector(0, z);
} else {
p.combineToSelf(rhoNext / rhoPrev, 1., z);
}
q = a.operate(p);
final double pq = p.dotProduct(q);
if (check && (pq <= 0.)) {
final NonPositiveDefiniteOperatorException e;
e = new NonPositiveDefiniteOperatorException();
final ExceptionContext context = e.getContext();
context.setValue(OPERATOR, a);
context.setValue(VECTOR, p);
throw e;
}
final double alpha = rhoNext / pq;
x.combineToSelf(1., alpha, p);
r.combineToSelf(1., -alpha, q);
rhoPrev = rhoNext;
rnorm = r.getNorm();
evt = new DefaultIterativeLinearSolverEvent(this,
manager.getIterations(), xro, bro, rro, rnorm);
manager.fireIterationPerformedEvent(evt);
if (rnorm <= rmax) {
manager.fireTerminationEvent(evt);
return x;
}
}
}
}