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 }