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 */ 017package org.apache.commons.math4.legacy.linear; 018 019import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 020import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 021import org.apache.commons.math4.legacy.exception.NullArgumentException; 022import org.apache.commons.math4.legacy.exception.util.ExceptionContext; 023 024/** 025 * <p> 026 * This is an implementation of the conjugate gradient method for 027 * {@link RealLinearOperator}. It follows closely the template by <a 028 * href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at 029 * hand is A · x = b, and the residual is r = b - A · x. 030 * </p> 031 * <p><b><a id="stopcrit">Default stopping criterion</a></b></p> 032 * <p> 033 * A default stopping criterion is implemented. The iterations stop when || r || 034 * ≤ δ || b ||, where b is the right-hand side vector, r the current 035 * estimate of the residual, and δ a user-specified tolerance. It should 036 * be noted that r is the so-called <em>updated</em> residual, which might 037 * differ from the true residual due to rounding-off errors (see e.g. <a 038 * href="#STRA2002">Strakos and Tichy, 2002</a>). 039 * </p> 040 * <p><b>Iteration count</b></p> 041 * <p> 042 * In the present context, an iteration should be understood as one evaluation 043 * of the matrix-vector product A · x. The initialization phase therefore 044 * counts as one iteration. 045 * </p> 046 * <p><b><a id="context">Exception context</a></b></p> 047 * <p> 048 * Besides standard {@link DimensionMismatchException}, this class might throw 049 * {@link NonPositiveDefiniteOperatorException} if the linear operator or 050 * the preconditioner are not positive definite. In this case, the 051 * {@link ExceptionContext} provides some more information 052 * <ul> 053 * <li>key {@code "operator"} points to the offending linear operator, say L,</li> 054 * <li>key {@code "vector"} points to the offending vector, say x, such that 055 * x<sup>T</sup> · L · x < 0.</li> 056 * </ul> 057 * 058 * <p><b>References</b></p> 059 * <dl> 060 * <dt><a id="BARR1994">Barret et al. (1994)</a></dt> 061 * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, 062 * V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst, 063 * <a href="http://www.netlib.org/linalg/html_templates/Templates.html"><em> 064 * Templates for the Solution of Linear Systems: Building Blocks for Iterative 065 * Methods</em></a>, SIAM</dd> 066 * <dt><a id="STRA2002">Strakos and Tichy (2002)</a></dt> 067 * <dd>Z. Strakos and P. Tichy, <a 068 * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf"> 069 * <em>On error estimation in the conjugate gradient method and why it works 070 * in finite precision computations</em></a>, Electronic Transactions on 071 * Numerical Analysis 13: 56-80, 2002</dd> 072 * </dl> 073 * 074 * @since 3.0 075 */ 076public class ConjugateGradient 077 extends PreconditionedIterativeLinearSolver { 078 079 /** Key for the <a href="#context">exception context</a>. */ 080 public static final String OPERATOR = "operator"; 081 082 /** Key for the <a href="#context">exception context</a>. */ 083 public static final String VECTOR = "vector"; 084 085 /** 086 * {@code true} if positive-definiteness of matrix and preconditioner should 087 * be checked. 088 */ 089 private boolean check; 090 091 /** The value of δ, for the default stopping criterion. */ 092 private final double delta; 093 094 /** 095 * Creates a new instance of this class, with <a href="#stopcrit">default 096 * stopping criterion</a>. 097 * 098 * @param maxIterations the maximum number of iterations 099 * @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}