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 agents.anac.y2019.harddealer.math3.linear;
|
---|
18 |
|
---|
19 | import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
|
---|
20 | import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
|
---|
21 | import agents.anac.y2019.harddealer.math3.exception.NullArgumentException;
|
---|
22 | import agents.anac.y2019.harddealer.math3.exception.util.ExceptionContext;
|
---|
23 | import agents.anac.y2019.harddealer.math3.util.IterationManager;
|
---|
24 |
|
---|
25 | /**
|
---|
26 | * <p>
|
---|
27 | * This is an implementation of the conjugate gradient method for
|
---|
28 | * {@link RealLinearOperator}. It follows closely the template by <a
|
---|
29 | * href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at
|
---|
30 | * hand is A · x = b, and the residual is r = b - A · x.
|
---|
31 | * </p>
|
---|
32 | * <h3><a id="stopcrit">Default stopping criterion</a></h3>
|
---|
33 | * <p>
|
---|
34 | * A default stopping criterion is implemented. The iterations stop when || r ||
|
---|
35 | * ≤ δ || b ||, where b is the right-hand side vector, r the current
|
---|
36 | * estimate of the residual, and δ a user-specified tolerance. It should
|
---|
37 | * be noted that r is the so-called <em>updated</em> residual, which might
|
---|
38 | * differ from the true residual due to rounding-off errors (see e.g. <a
|
---|
39 | * href="#STRA2002">Strakos and Tichy, 2002</a>).
|
---|
40 | * </p>
|
---|
41 | * <h3>Iteration count</h3>
|
---|
42 | * <p>
|
---|
43 | * In the present context, an iteration should be understood as one evaluation
|
---|
44 | * of the matrix-vector product A · x. The initialization phase therefore
|
---|
45 | * counts as one iteration.
|
---|
46 | * </p>
|
---|
47 | * <h3><a id="context">Exception context</a></h3>
|
---|
48 | * <p>
|
---|
49 | * Besides standard {@link DimensionMismatchException}, this class might throw
|
---|
50 | * {@link NonPositiveDefiniteOperatorException} if the linear operator or
|
---|
51 | * the preconditioner are not positive definite. In this case, the
|
---|
52 | * {@link ExceptionContext} provides some more information
|
---|
53 | * <ul>
|
---|
54 | * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
|
---|
55 | * <li>key {@code "vector"} points to the offending vector, say x, such that
|
---|
56 | * x<sup>T</sup> · L · x < 0.</li>
|
---|
57 | * </ul>
|
---|
58 | * </p>
|
---|
59 | * <h3>References</h3>
|
---|
60 | * <dl>
|
---|
61 | * <dt><a id="BARR1994">Barret et al. (1994)</a></dt>
|
---|
62 | * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra,
|
---|
63 | * V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
|
---|
64 | * <a href="http://www.netlib.org/linalg/html_templates/Templates.html"><em>
|
---|
65 | * Templates for the Solution of Linear Systems: Building Blocks for Iterative
|
---|
66 | * Methods</em></a>, SIAM</dd>
|
---|
67 | * <dt><a id="STRA2002">Strakos and Tichy (2002)
|
---|
68 | * <dt>
|
---|
69 | * <dd>Z. Strakos and P. Tichy, <a
|
---|
70 | * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf">
|
---|
71 | * <em>On error estimation in the conjugate gradient method and why it works
|
---|
72 | * in finite precision computations</em></a>, Electronic Transactions on
|
---|
73 | * Numerical Analysis 13: 56-80, 2002</dd>
|
---|
74 | * </dl>
|
---|
75 | *
|
---|
76 | * @since 3.0
|
---|
77 | */
|
---|
78 | public class ConjugateGradient
|
---|
79 | extends PreconditionedIterativeLinearSolver {
|
---|
80 |
|
---|
81 | /** Key for the <a href="#context">exception context</a>. */
|
---|
82 | public static final String OPERATOR = "operator";
|
---|
83 |
|
---|
84 | /** Key for the <a href="#context">exception context</a>. */
|
---|
85 | public static final String VECTOR = "vector";
|
---|
86 |
|
---|
87 | /**
|
---|
88 | * {@code true} if positive-definiteness of matrix and preconditioner should
|
---|
89 | * be checked.
|
---|
90 | */
|
---|
91 | private boolean check;
|
---|
92 |
|
---|
93 | /** The value of δ, for the default stopping criterion. */
|
---|
94 | private final double delta;
|
---|
95 |
|
---|
96 | /**
|
---|
97 | * Creates a new instance of this class, with <a href="#stopcrit">default
|
---|
98 | * stopping criterion</a>.
|
---|
99 | *
|
---|
100 | * @param maxIterations the maximum number of iterations
|
---|
101 | * @param delta the δ parameter for the default stopping criterion
|
---|
102 | * @param check {@code true} if positive definiteness of both matrix and
|
---|
103 | * preconditioner should be checked
|
---|
104 | */
|
---|
105 | public ConjugateGradient(final int maxIterations, final double delta,
|
---|
106 | final boolean check) {
|
---|
107 | super(maxIterations);
|
---|
108 | this.delta = delta;
|
---|
109 | this.check = check;
|
---|
110 | }
|
---|
111 |
|
---|
112 | /**
|
---|
113 | * Creates a new instance of this class, with <a href="#stopcrit">default
|
---|
114 | * stopping criterion</a> and custom iteration manager.
|
---|
115 | *
|
---|
116 | * @param manager the custom iteration manager
|
---|
117 | * @param delta the δ parameter for the default stopping criterion
|
---|
118 | * @param check {@code true} if positive definiteness of both matrix and
|
---|
119 | * preconditioner should be checked
|
---|
120 | * @throws NullArgumentException if {@code manager} is {@code null}
|
---|
121 | */
|
---|
122 | public ConjugateGradient(final IterationManager manager,
|
---|
123 | final double delta, final boolean check)
|
---|
124 | throws NullArgumentException {
|
---|
125 | super(manager);
|
---|
126 | this.delta = delta;
|
---|
127 | this.check = check;
|
---|
128 | }
|
---|
129 |
|
---|
130 | /**
|
---|
131 | * Returns {@code true} if positive-definiteness should be checked for both
|
---|
132 | * matrix and preconditioner.
|
---|
133 | *
|
---|
134 | * @return {@code true} if the tests are to be performed
|
---|
135 | */
|
---|
136 | public final boolean getCheck() {
|
---|
137 | return check;
|
---|
138 | }
|
---|
139 |
|
---|
140 | /**
|
---|
141 | * {@inheritDoc}
|
---|
142 | *
|
---|
143 | * @throws NonPositiveDefiniteOperatorException if {@code a} or {@code m} is
|
---|
144 | * not positive definite
|
---|
145 | */
|
---|
146 | @Override
|
---|
147 | public RealVector solveInPlace(final RealLinearOperator a,
|
---|
148 | final RealLinearOperator m,
|
---|
149 | final RealVector b,
|
---|
150 | final RealVector x0)
|
---|
151 | throws NullArgumentException, NonPositiveDefiniteOperatorException,
|
---|
152 | NonSquareOperatorException, DimensionMismatchException,
|
---|
153 | MaxCountExceededException {
|
---|
154 | checkParameters(a, m, b, x0);
|
---|
155 | final IterationManager manager = getIterationManager();
|
---|
156 | // Initialization of default stopping criterion
|
---|
157 | manager.resetIterationCount();
|
---|
158 | final double rmax = delta * b.getNorm();
|
---|
159 | final RealVector bro = RealVector.unmodifiableRealVector(b);
|
---|
160 |
|
---|
161 | // Initialization phase counts as one iteration.
|
---|
162 | manager.incrementIterationCount();
|
---|
163 | // p and x are constructed as copies of x0, since presumably, the type
|
---|
164 | // of x is optimized for the calculation of the matrix-vector product
|
---|
165 | // A.x.
|
---|
166 | final RealVector x = x0;
|
---|
167 | final RealVector xro = RealVector.unmodifiableRealVector(x);
|
---|
168 | final RealVector p = x.copy();
|
---|
169 | RealVector q = a.operate(p);
|
---|
170 |
|
---|
171 | final RealVector r = b.combine(1, -1, q);
|
---|
172 | final RealVector rro = RealVector.unmodifiableRealVector(r);
|
---|
173 | double rnorm = r.getNorm();
|
---|
174 | RealVector z;
|
---|
175 | if (m == null) {
|
---|
176 | z = r;
|
---|
177 | } else {
|
---|
178 | z = null;
|
---|
179 | }
|
---|
180 | IterativeLinearSolverEvent evt;
|
---|
181 | evt = new DefaultIterativeLinearSolverEvent(this,
|
---|
182 | manager.getIterations(), xro, bro, rro, rnorm);
|
---|
183 | manager.fireInitializationEvent(evt);
|
---|
184 | if (rnorm <= rmax) {
|
---|
185 | manager.fireTerminationEvent(evt);
|
---|
186 | return x;
|
---|
187 | }
|
---|
188 | double rhoPrev = 0.;
|
---|
189 | while (true) {
|
---|
190 | manager.incrementIterationCount();
|
---|
191 | evt = new DefaultIterativeLinearSolverEvent(this,
|
---|
192 | manager.getIterations(), xro, bro, rro, rnorm);
|
---|
193 | manager.fireIterationStartedEvent(evt);
|
---|
194 | if (m != null) {
|
---|
195 | z = m.operate(r);
|
---|
196 | }
|
---|
197 | final double rhoNext = r.dotProduct(z);
|
---|
198 | if (check && (rhoNext <= 0.)) {
|
---|
199 | final NonPositiveDefiniteOperatorException e;
|
---|
200 | e = new NonPositiveDefiniteOperatorException();
|
---|
201 | final ExceptionContext context = e.getContext();
|
---|
202 | context.setValue(OPERATOR, m);
|
---|
203 | context.setValue(VECTOR, r);
|
---|
204 | throw e;
|
---|
205 | }
|
---|
206 | if (manager.getIterations() == 2) {
|
---|
207 | p.setSubVector(0, z);
|
---|
208 | } else {
|
---|
209 | p.combineToSelf(rhoNext / rhoPrev, 1., z);
|
---|
210 | }
|
---|
211 | q = a.operate(p);
|
---|
212 | final double pq = p.dotProduct(q);
|
---|
213 | if (check && (pq <= 0.)) {
|
---|
214 | final NonPositiveDefiniteOperatorException e;
|
---|
215 | e = new NonPositiveDefiniteOperatorException();
|
---|
216 | final ExceptionContext context = e.getContext();
|
---|
217 | context.setValue(OPERATOR, a);
|
---|
218 | context.setValue(VECTOR, p);
|
---|
219 | throw e;
|
---|
220 | }
|
---|
221 | final double alpha = rhoNext / pq;
|
---|
222 | x.combineToSelf(1., alpha, p);
|
---|
223 | r.combineToSelf(1., -alpha, q);
|
---|
224 | rhoPrev = rhoNext;
|
---|
225 | rnorm = r.getNorm();
|
---|
226 | evt = new DefaultIterativeLinearSolverEvent(this,
|
---|
227 | manager.getIterations(), xro, bro, rro, rnorm);
|
---|
228 | manager.fireIterationPerformedEvent(evt);
|
---|
229 | if (rnorm <= rmax) {
|
---|
230 | manager.fireTerminationEvent(evt);
|
---|
231 | return x;
|
---|
232 | }
|
---|
233 | }
|
---|
234 | }
|
---|
235 | }
|
---|