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; 18 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; 23 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 · x = b, and the residual is r = b - A · 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 * ≤ δ || b ||, where b is the right-hand side vector, r the current 35 * estimate of the residual, and δ 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 · 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> · L · x < 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 { 78 79 /** Key for the <a href="#context">exception context</a>. */ 80 public static final String OPERATOR = "operator"; 81 82 /** Key for the <a href="#context">exception context</a>. */ 83 public static final String VECTOR = "vector"; 84 85 /** 86 * {@code true} if positive-definiteness of matrix and preconditioner should 87 * be checked. 88 */ 89 private boolean check; 90 91 /** The value of δ, for the default stopping criterion. */ 92 private final double delta; 93 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 δ 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 } 109 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 δ 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 } 127 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 } 137 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); 158 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); 168 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 }