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.FastMath;
024    import org.apache.commons.math3.util.IterationManager;
025    import org.apache.commons.math3.util.MathUtils;
026    
027    /**
028     * <p>
029     * Implementation of the SYMMLQ iterative linear solver proposed by <a
030     * href="#PAIG1975">Paige and Saunders (1975)</a>. This implementation is
031     * largely based on the FORTRAN code by Pr. Michael A. Saunders, available <a
032     * href="http://www.stanford.edu/group/SOL/software/symmlq/f77/">here</a>.
033     * </p>
034     * <p>
035     * SYMMLQ is designed to solve the system of linear equations A &middot; x = b
036     * where A is an n &times; n self-adjoint linear operator (defined as a
037     * {@link RealLinearOperator}), and b is a given vector. The operator A is not
038     * required to be positive definite. If A is known to be definite, the method of
039     * conjugate gradients might be preferred, since it will require about the same
040     * number of iterations as SYMMLQ but slightly less work per iteration.
041     * </p>
042     * <p>
043     * SYMMLQ is designed to solve the system (A - shift &middot; I) &middot; x = b,
044     * where shift is a specified scalar value. If shift and b are suitably chosen,
045     * the computed vector x may approximate an (unnormalized) eigenvector of A, as
046     * in the methods of inverse iteration and/or Rayleigh-quotient iteration.
047     * Again, the linear operator (A - shift &middot; I) need not be positive
048     * definite (but <em>must</em> be self-adjoint). The work per iteration is very
049     * slightly less if shift = 0.
050     * </p>
051     * <h3>Preconditioning</h3>
052     * <p>
053     * Preconditioning may reduce the number of iterations required. The solver may
054     * be provided with a positive definite preconditioner
055     * M = P<sup>T</sup> &middot; P
056     * that is known to approximate
057     * (A - shift &middot; I)<sup>-1</sup> in some sense, where matrix-vector
058     * products of the form M &middot; y = x can be computed efficiently. Then
059     * SYMMLQ will implicitly solve the system of equations
060     * P &middot; (A - shift &middot; I) &middot; P<sup>T</sup> &middot;
061     * x<sub>hat</sub> = P &middot; b, i.e.
062     * A<sub>hat</sub> &middot; x<sub>hat</sub> = b<sub>hat</sub>,
063     * where
064     * A<sub>hat</sub> = P &middot; (A - shift &middot; I) &middot; P<sup>T</sup>,
065     * b<sub>hat</sub> = P &middot; b,
066     * and return the solution
067     * x = P<sup>T</sup> &middot; x<sub>hat</sub>.
068     * The associated residual is
069     * r<sub>hat</sub> = b<sub>hat</sub> - A<sub>hat</sub> &middot; x<sub>hat</sub>
070     *                 = P &middot; [b - (A - shift &middot; I) &middot; x]
071     *                 = P &middot; r.
072     * </p>
073     * <p>
074     * In the case of preconditioning, the {@link IterativeLinearSolverEvent}s that
075     * this solver fires are such that
076     * {@link IterativeLinearSolverEvent#getNormOfResidual()} returns the norm of
077     * the <em>preconditioned</em>, updated residual, ||P &middot; r||, not the norm
078     * of the <em>true</em> residual ||r||.
079     * </p>
080     * <h3><a id="stopcrit">Default stopping criterion</a></h3>
081     * <p>
082     * A default stopping criterion is implemented. The iterations stop when || rhat
083     * || &le; &delta; || Ahat || || xhat ||, where xhat is the current estimate of
084     * the solution of the transformed system, rhat the current estimate of the
085     * corresponding residual, and &delta; a user-specified tolerance.
086     * </p>
087     * <h3>Iteration count</h3>
088     * <p>
089     * In the present context, an iteration should be understood as one evaluation
090     * of the matrix-vector product A &middot; x. The initialization phase therefore
091     * counts as one iteration. If the user requires checks on the symmetry of A,
092     * this entails one further matrix-vector product in the initial phase. This
093     * further product is <em>not</em> accounted for in the iteration count. In
094     * other words, the number of iterations required to reach convergence will be
095     * identical, whether checks have been required or not.
096     * </p>
097     * <p>
098     * The present definition of the iteration count differs from that adopted in
099     * the original FOTRAN code, where the initialization phase was <em>not</em>
100     * taken into account.
101     * </p>
102     * <h3><a id="initguess">Initial guess of the solution</a></h3>
103     * <p>
104     * The {@code x} parameter in
105     * <ul>
106     * <li>{@link #solve(RealLinearOperator, RealVector, RealVector)},</li>
107     * <li>{@link #solve(RealLinearOperator, RealLinearOperator, RealVector, RealVector)}},</li>
108     * <li>{@link #solveInPlace(RealLinearOperator, RealVector, RealVector)},</li>
109     * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector)},</li>
110     * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector, boolean, double)},</li>
111     * </ul>
112     * should not be considered as an initial guess, as it is set to zero in the
113     * initial phase. If x<sub>0</sub> is known to be a good approximation to x, one
114     * should compute r<sub>0</sub> = b - A &middot; x, solve A &middot; dx = r0,
115     * and set x = x<sub>0</sub> + dx.
116     * </p>
117     * <h3><a id="context">Exception context</a></h3>
118     * <p>
119     * Besides standard {@link DimensionMismatchException}, this class might throw
120     * {@link NonSelfAdjointOperatorException} if the linear operator or the
121     * preconditioner are not symmetric. In this case, the {@link ExceptionContext}
122     * provides more information
123     * <ul>
124     * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
125     * <li>key {@code "vector1"} points to the first offending vector, say x,
126     * <li>key {@code "vector2"} points to the second offending vector, say y, such
127     * that x<sup>T</sup> &middot; L &middot; y &ne; y<sup>T</sup> &middot; L
128     * &middot; x (within a certain accuracy).</li>
129     * </ul>
130     * </p>
131     * <p>
132     * {@link NonPositiveDefiniteOperatorException} might also be thrown in case the
133     * preconditioner is not positive definite. The relevant keys to the
134     * {@link ExceptionContext} are
135     * <ul>
136     * <li>key {@code "operator"}, which points to the offending linear operator,
137     * say L,</li>
138     * <li>key {@code "vector"}, which points to the offending vector, say x, such
139     * that x<sup>T</sup> &middot; L &middot; x < 0.</li>
140     * </ul>
141     * </p>
142     * <h3>References</h3>
143     * <dl>
144     * <dt><a id="PAIG1975">Paige and Saunders (1975)</a></dt>
145     * <dd>C. C. Paige and M. A. Saunders, <a
146     * href="http://www.stanford.edu/group/SOL/software/symmlq/PS75.pdf"><em>
147     * Solution of Sparse Indefinite Systems of Linear Equations</em></a>, SIAM
148     * Journal on Numerical Analysis 12(4): 617-629, 1975</dd>
149     * </dl>
150     *
151     * @version $Id: SymmLQ.java 1416643 2012-12-03 19:37:14Z tn $
152     * @since 3.0
153     */
154    public class SymmLQ
155        extends PreconditionedIterativeLinearSolver {
156    
157        /*
158         * IMPLEMENTATION NOTES
159         * --------------------
160         * The implementation follows as closely as possible the notations of Paige
161         * and Saunders (1975). Attention must be paid to the fact that some
162         * quantities which are relevant to iteration k can only be computed in
163         * iteration (k+1). Therefore, minute attention must be paid to the index of
164         * each state variable of this algorithm.
165         *
166         * 1. Preconditioning
167         *    ---------------
168         * The Lanczos iterations associated with Ahat and bhat read
169         *   beta[1] = ||P * b||
170         *   v[1] = P * b / beta[1]
171         *   beta[k+1] * v[k+1] = Ahat * v[k] - alpha[k] * v[k] - beta[k] * v[k-1]
172         *                      = P * (A - shift * I) * P' * v[k] - alpha[k] * v[k]
173         *                        - beta[k] * v[k-1]
174         * Multiplying both sides by P', we get
175         *   beta[k+1] * (P' * v)[k+1] = M * (A - shift * I) * (P' * v)[k]
176         *                               - alpha[k] * (P' * v)[k]
177         *                               - beta[k] * (P' * v[k-1]),
178         * and
179         *   alpha[k+1] = v[k+1]' * Ahat * v[k+1]
180         *              = v[k+1]' * P * (A - shift * I) * P' * v[k+1]
181         *              = (P' * v)[k+1]' * (A - shift * I) * (P' * v)[k+1].
182         *
183         * In other words, the Lanczos iterations are unchanged, except for the fact
184         * that we really compute (P' * v) instead of v. It can easily be checked
185         * that all other formulas are unchanged. It must be noted that P is never
186         * explicitly used, only matrix-vector products involving are invoked.
187         *
188         * 2. Accounting for the shift parameter
189         *    ----------------------------------
190         * Is trivial: each time A.operate(x) is invoked, one must subtract shift * x
191         * to the result.
192         *
193         * 3. Accounting for the goodb flag
194         *    -----------------------------
195         * When goodb is set to true, the component of xL along b is computed
196         * separately. From Paige and Saunders (1975), equation (5.9), we have
197         *   wbar[k+1] = s[k] * wbar[k] - c[k] * v[k+1],
198         *   wbar[1] = v[1].
199         * Introducing wbar2[k] = wbar[k] - s[1] * ... * s[k-1] * v[1], it can
200         * easily be verified by induction that wbar2 follows the same recursive
201         * relation
202         *   wbar2[k+1] = s[k] * wbar2[k] - c[k] * v[k+1],
203         *   wbar2[1] = 0,
204         * and we then have
205         *   w[k] = c[k] * wbar2[k] + s[k] * v[k+1]
206         *          + s[1] * ... * s[k-1] * c[k] * v[1].
207         * Introducing w2[k] = w[k] - s[1] * ... * s[k-1] * c[k] * v[1], we find,
208         * from (5.10)
209         *   xL[k] = zeta[1] * w[1] + ... + zeta[k] * w[k]
210         *         = zeta[1] * w2[1] + ... + zeta[k] * w2[k]
211         *           + (s[1] * c[2] * zeta[2] + ...
212         *           + s[1] * ... * s[k-1] * c[k] * zeta[k]) * v[1]
213         *         = xL2[k] + bstep[k] * v[1],
214         * where xL2[k] is defined by
215         *   xL2[0] = 0,
216         *   xL2[k+1] = xL2[k] + zeta[k+1] * w2[k+1],
217         * and bstep is defined by
218         *   bstep[1] = 0,
219         *   bstep[k] = bstep[k-1] + s[1] * ... * s[k-1] * c[k] * zeta[k].
220         * We also have, from (5.11)
221         *   xC[k] = xL[k-1] + zbar[k] * wbar[k]
222         *         = xL2[k-1] + zbar[k] * wbar2[k]
223         *           + (bstep[k-1] + s[1] * ... * s[k-1] * zbar[k]) * v[1].
224         */
225    
226        /**
227         * <p>
228         * A simple container holding the non-final variables used in the
229         * iterations. Making the current state of the solver visible from the
230         * outside is necessary, because during the iterations, {@code x} does not
231         * <em>exactly</em> hold the current estimate of the solution. Indeed,
232         * {@code x} needs in general to be moved from the LQ point to the CG point.
233         * Besides, additional upudates must be carried out in case {@code goodb} is
234         * set to {@code true}.
235         * </p>
236         * <p>
237         * In all subsequent comments, the description of the state variables refer
238         * to their value after a call to {@link #update()}. In these comments, k is
239         * the current number of evaluations of matrix-vector products.
240         * </p>
241         */
242        private static class State {
243            /** The cubic root of {@link #MACH_PREC}. */
244            static final double CBRT_MACH_PREC;
245    
246            /** The machine precision. */
247            static final double MACH_PREC;
248    
249            /** Reference to the linear operator. */
250            private final RealLinearOperator a;
251    
252            /** Reference to the right-hand side vector. */
253            private final RealVector b;
254    
255            /** {@code true} if symmetry of matrix and conditioner must be checked. */
256            private final boolean check;
257    
258            /**
259             * The value of the custom tolerance &delta; for the default stopping
260             * criterion.
261             */
262            private final double delta;
263    
264            /** The value of beta[k+1]. */
265            private double beta;
266    
267            /** The value of beta[1]. */
268            private double beta1;
269    
270            /** The value of bstep[k-1]. */
271            private double bstep;
272    
273            /** The estimate of the norm of P * rC[k]. */
274            private double cgnorm;
275    
276            /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */
277            private double dbar;
278    
279            /**
280             * The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the
281             * initial code.
282             */
283            private double gammaZeta;
284    
285            /** The value of gbar[k]. */
286            private double gbar;
287    
288            /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
289            private double gmax;
290    
291            /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
292            private double gmin;
293    
294            /** Copy of the {@code goodb} parameter. */
295            private final boolean goodb;
296    
297            /** {@code true} if the default convergence criterion is verified. */
298            private boolean hasConverged;
299    
300            /** The estimate of the norm of P * rL[k-1]. */
301            private double lqnorm;
302    
303            /** Reference to the preconditioner, M. */
304            private final RealLinearOperator m;
305    
306            /**
307             * The value of (-eps[k+1] * zeta[k-1]). Was called {@code rhs2} in the
308             * initial code.
309             */
310            private double minusEpsZeta;
311    
312            /** The value of M * b. */
313            private final RealVector mb;
314    
315            /** The value of beta[k]. */
316            private double oldb;
317    
318            /** The value of beta[k] * M^(-1) * P' * v[k]. */
319            private RealVector r1;
320    
321            /** The value of beta[k+1] * M^(-1) * P' * v[k+1]. */
322            private RealVector r2;
323    
324            /**
325             * The value of the updated, preconditioned residual P * r. This value is
326             * given by {@code min(}{@link #cgnorm}{@code , }{@link #lqnorm}{@code )}.
327             */
328            private double rnorm;
329    
330            /** Copy of the {@code shift} parameter. */
331            private final double shift;
332    
333            /** The value of s[1] * ... * s[k-1]. */
334            private double snprod;
335    
336            /**
337             * An estimate of the square of the norm of A * V[k], based on Paige and
338             * Saunders (1975), equation (3.3).
339             */
340            private double tnorm;
341    
342            /**
343             * The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] *
344             * v[1]) if {@code goodb} is {@code true}. Was called {@code w} in the
345             * initial code.
346             */
347            private RealVector wbar;
348    
349            /**
350             * A reference to the vector to be updated with the solution. Contains
351             * the value of xL[k-1] if {@code goodb} is {@code false}, (xL[k-1] -
352             * bstep[k-1] * v[1]) otherwise.
353             */
354            private final RealVector xL;
355    
356            /** The value of beta[k+1] * P' * v[k+1]. */
357            private RealVector y;
358    
359            /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */
360            private double ynorm2;
361    
362            /** The value of {@code b == 0} (exact floating-point equality). */
363            private boolean bIsNull;
364    
365            static {
366                MACH_PREC = FastMath.ulp(1.);
367                CBRT_MACH_PREC = FastMath.cbrt(MACH_PREC);
368            }
369    
370            /**
371             * Creates and inits to k = 1 a new instance of this class.
372             *
373             * @param a the linear operator A of the system
374             * @param m the preconditioner, M (can be {@code null})
375             * @param b the right-hand side vector
376             * @param goodb usually {@code false}, except if {@code x} is expected
377             * to contain a large multiple of {@code b}
378             * @param shift the amount to be subtracted to all diagonal elements of
379             * A
380             * @param delta the &delta; parameter for the default stopping criterion
381             * @param check {@code true} if self-adjointedness of both matrix and
382             * preconditioner should be checked
383             */
384            public State(final RealLinearOperator a,
385                final RealLinearOperator m,
386                final RealVector b,
387                final boolean goodb,
388                final double shift,
389                final double delta,
390                final boolean check) {
391                this.a = a;
392                this.m = m;
393                this.b = b;
394                this.xL = new ArrayRealVector(b.getDimension());
395                this.goodb = goodb;
396                this.shift = shift;
397                this.mb = m == null ? b : m.operate(b);
398                this.hasConverged = false;
399                this.check = check;
400                this.delta = delta;
401            }
402    
403            /**
404             * Performs a symmetry check on the specified linear operator, and throws an
405             * exception in case this check fails. Given a linear operator L, and a
406             * vector x, this method checks that
407             * x' &middot; L &middot; y = y' &middot; L &middot; x
408             * (within a given accuracy), where y = L &middot; x.
409             *
410             * @param l the linear operator L
411             * @param x the candidate vector x
412             * @param y the candidate vector y = L &middot; x
413             * @param z the vector z = L &middot; y
414             * @throws NonSelfAdjointOperatorException when the test fails
415             */
416            private static void checkSymmetry(final RealLinearOperator l,
417                final RealVector x, final RealVector y, final RealVector z)
418                throws NonSelfAdjointOperatorException {
419                final double s = y.dotProduct(y);
420                final double t = x.dotProduct(z);
421                final double epsa = (s + MACH_PREC) * CBRT_MACH_PREC;
422                if (FastMath.abs(s - t) > epsa) {
423                    final NonSelfAdjointOperatorException e;
424                    e = new NonSelfAdjointOperatorException();
425                    final ExceptionContext context = e.getContext();
426                    context.setValue(SymmLQ.OPERATOR, l);
427                    context.setValue(SymmLQ.VECTOR1, x);
428                    context.setValue(SymmLQ.VECTOR2, y);
429                    context.setValue(SymmLQ.THRESHOLD, Double.valueOf(epsa));
430                    throw e;
431                }
432            }
433    
434            /**
435             * Throws a new {@link NonPositiveDefiniteOperatorException} with
436             * appropriate context.
437             *
438             * @param l the offending linear operator
439             * @param v the offending vector
440             * @throws NonPositiveDefiniteOperatorException in any circumstances
441             */
442            private static void throwNPDLOException(final RealLinearOperator l,
443                final RealVector v) throws NonPositiveDefiniteOperatorException {
444                final NonPositiveDefiniteOperatorException e;
445                e = new NonPositiveDefiniteOperatorException();
446                final ExceptionContext context = e.getContext();
447                context.setValue(OPERATOR, l);
448                context.setValue(VECTOR, v);
449                throw e;
450            }
451    
452            /**
453             * A clone of the BLAS {@code DAXPY} function, which carries out the
454             * operation y &larr; a &middot; x + y. This is for internal use only: no
455             * dimension checks are provided.
456             *
457             * @param a the scalar by which {@code x} is to be multiplied
458             * @param x the vector to be added to {@code y}
459             * @param y the vector to be incremented
460             */
461            private static void daxpy(final double a, final RealVector x,
462                final RealVector y) {
463                final int n = x.getDimension();
464                for (int i = 0; i < n; i++) {
465                    y.setEntry(i, a * x.getEntry(i) + y.getEntry(i));
466                }
467            }
468    
469            /**
470             * A BLAS-like function, for the operation z &larr; a &middot; x + b
471             * &middot; y + z. This is for internal use only: no dimension checks are
472             * provided.
473             *
474             * @param a the scalar by which {@code x} is to be multiplied
475             * @param x the first vector to be added to {@code z}
476             * @param b the scalar by which {@code y} is to be multiplied
477             * @param y the second vector to be added to {@code z}
478             * @param z the vector to be incremented
479             */
480            private static void daxpbypz(final double a, final RealVector x,
481                final double b, final RealVector y, final RealVector z) {
482                final int n = z.getDimension();
483                for (int i = 0; i < n; i++) {
484                    final double zi;
485                    zi = a * x.getEntry(i) + b * y.getEntry(i) + z.getEntry(i);
486                    z.setEntry(i, zi);
487                }
488            }
489    
490            /**
491             * <p>
492             * Move to the CG point if it seems better. In this version of SYMMLQ,
493             * the convergence tests involve only cgnorm, so we're unlikely to stop
494             * at an LQ point, except if the iteration limit interferes.
495             * </p>
496             * <p>
497             * Additional upudates are also carried out in case {@code goodb} is set
498             * to {@code true}.
499             * </p>
500             *
501             * @param x the vector to be updated with the refined value of xL
502             */
503             void refineSolution(final RealVector x) {
504                final int n = this.xL.getDimension();
505                if (lqnorm < cgnorm) {
506                    if (!goodb) {
507                        x.setSubVector(0, this.xL);
508                    } else {
509                        final double step = bstep / beta1;
510                        for (int i = 0; i < n; i++) {
511                            final double bi = mb.getEntry(i);
512                            final double xi = this.xL.getEntry(i);
513                            x.setEntry(i, xi + step * bi);
514                        }
515                    }
516                } else {
517                    final double anorm = FastMath.sqrt(tnorm);
518                    final double diag = gbar == 0. ? anorm * MACH_PREC : gbar;
519                    final double zbar = gammaZeta / diag;
520                    final double step = (bstep + snprod * zbar) / beta1;
521                    // ynorm = FastMath.sqrt(ynorm2 + zbar * zbar);
522                    if (!goodb) {
523                        for (int i = 0; i < n; i++) {
524                            final double xi = this.xL.getEntry(i);
525                            final double wi = wbar.getEntry(i);
526                            x.setEntry(i, xi + zbar * wi);
527                        }
528                    } else {
529                        for (int i = 0; i < n; i++) {
530                            final double xi = this.xL.getEntry(i);
531                            final double wi = wbar.getEntry(i);
532                            final double bi = mb.getEntry(i);
533                            x.setEntry(i, xi + zbar * wi + step * bi);
534                        }
535                    }
536                }
537            }
538    
539            /**
540             * Performs the initial phase of the SYMMLQ algorithm. On return, the
541             * value of the state variables of {@code this} object correspond to k =
542             * 1.
543             */
544             void init() {
545                this.xL.set(0.);
546                /*
547                 * Set up y for the first Lanczos vector. y and beta1 will be zero
548                 * if b = 0.
549                 */
550                this.r1 = this.b.copy();
551                this.y = this.m == null ? this.b.copy() : this.m.operate(this.r1);
552                if ((this.m != null) && this.check) {
553                    checkSymmetry(this.m, this.r1, this.y, this.m.operate(this.y));
554                }
555    
556                this.beta1 = this.r1.dotProduct(this.y);
557                if (this.beta1 < 0.) {
558                    throwNPDLOException(this.m, this.y);
559                }
560                if (this.beta1 == 0.) {
561                    /* If b = 0 exactly, stop with x = 0. */
562                    this.bIsNull = true;
563                    return;
564                }
565                this.bIsNull = false;
566                this.beta1 = FastMath.sqrt(this.beta1);
567                /* At this point
568                 *   r1 = b,
569                 *   y = M * b,
570                 *   beta1 = beta[1].
571                 */
572                final RealVector v = this.y.mapMultiply(1. / this.beta1);
573                this.y = this.a.operate(v);
574                if (this.check) {
575                    checkSymmetry(this.a, v, this.y, this.a.operate(this.y));
576                }
577                /*
578                 * Set up y for the second Lanczos vector. y and beta will be zero
579                 * or very small if b is an eigenvector.
580                 */
581                daxpy(-this.shift, v, this.y);
582                final double alpha = v.dotProduct(this.y);
583                daxpy(-alpha / this.beta1, this.r1, this.y);
584                /*
585                 * At this point
586                 *   alpha = alpha[1]
587                 *   y     = beta[2] * M^(-1) * P' * v[2]
588                 */
589                /* Make sure r2 will be orthogonal to the first v. */
590                final double vty = v.dotProduct(this.y);
591                final double vtv = v.dotProduct(v);
592                daxpy(-vty / vtv, v, this.y);
593                this.r2 = this.y.copy();
594                if (this.m != null) {
595                    this.y = this.m.operate(this.r2);
596                }
597                this.oldb = this.beta1;
598                this.beta = this.r2.dotProduct(this.y);
599                if (this.beta < 0.) {
600                    throwNPDLOException(this.m, this.y);
601                }
602                this.beta = FastMath.sqrt(this.beta);
603                /*
604                 * At this point
605                 *   oldb = beta[1]
606                 *   beta = beta[2]
607                 *   y  = beta[2] * P' * v[2]
608                 *   r2 = beta[2] * M^(-1) * P' * v[2]
609                 */
610                this.cgnorm = this.beta1;
611                this.gbar = alpha;
612                this.dbar = this.beta;
613                this.gammaZeta = this.beta1;
614                this.minusEpsZeta = 0.;
615                this.bstep = 0.;
616                this.snprod = 1.;
617                this.tnorm = alpha * alpha + this.beta * this.beta;
618                this.ynorm2 = 0.;
619                this.gmax = FastMath.abs(alpha) + MACH_PREC;
620                this.gmin = this.gmax;
621    
622                if (this.goodb) {
623                    this.wbar = new ArrayRealVector(this.a.getRowDimension());
624                    this.wbar.set(0.);
625                } else {
626                    this.wbar = v;
627                }
628                updateNorms();
629            }
630    
631            /**
632             * Performs the next iteration of the algorithm. The iteration count
633             * should be incremented prior to calling this method. On return, the
634             * value of the state variables of {@code this} object correspond to the
635             * current iteration count {@code k}.
636             */
637            void update() {
638                final RealVector v = y.mapMultiply(1. / beta);
639                y = a.operate(v);
640                daxpbypz(-shift, v, -beta / oldb, r1, y);
641                final double alpha = v.dotProduct(y);
642                /*
643                 * At this point
644                 *   v     = P' * v[k],
645                 *   y     = (A - shift * I) * P' * v[k] - beta[k] * M^(-1) * P' * v[k-1],
646                 *   alpha = v'[k] * P * (A - shift * I) * P' * v[k]
647                 *           - beta[k] * v[k]' * P * M^(-1) * P' * v[k-1]
648                 *         = v'[k] * P * (A - shift * I) * P' * v[k]
649                 *           - beta[k] * v[k]' * v[k-1]
650                 *         = alpha[k].
651                 */
652                daxpy(-alpha / beta, r2, y);
653                /*
654                 * At this point
655                 *   y = (A - shift * I) * P' * v[k] - alpha[k] * M^(-1) * P' * v[k]
656                 *       - beta[k] * M^(-1) * P' * v[k-1]
657                 *     = M^(-1) * P' * (P * (A - shift * I) * P' * v[k] -alpha[k] * v[k]
658                 *       - beta[k] * v[k-1])
659                 *     = beta[k+1] * M^(-1) * P' * v[k+1],
660                 * from Paige and Saunders (1975), equation (3.2).
661                 *
662                 * WATCH-IT: the two following lines work only because y is no longer
663                 * updated up to the end of the present iteration, and is
664                 * reinitialized at the beginning of the next iteration.
665                 */
666                r1 = r2;
667                r2 = y;
668                if (m != null) {
669                    y = m.operate(r2);
670                }
671                oldb = beta;
672                beta = r2.dotProduct(y);
673                if (beta < 0.) {
674                    throwNPDLOException(m, y);
675                }
676                beta = FastMath.sqrt(beta);
677                /*
678                 * At this point
679                 *   r1 = beta[k] * M^(-1) * P' * v[k],
680                 *   r2 = beta[k+1] * M^(-1) * P' * v[k+1],
681                 *   y  = beta[k+1] * P' * v[k+1],
682                 *   oldb = beta[k],
683                 *   beta = beta[k+1].
684                 */
685                tnorm += alpha * alpha + oldb * oldb + beta * beta;
686                /*
687                 * Compute the next plane rotation for Q. See Paige and Saunders
688                 * (1975), equation (5.6), with
689                 *   gamma = gamma[k-1],
690                 *   c     = c[k-1],
691                 *   s     = s[k-1].
692                 */
693                final double gamma = FastMath.sqrt(gbar * gbar + oldb * oldb);
694                final double c = gbar / gamma;
695                final double s = oldb / gamma;
696                /*
697                 * The relations
698                 *   gbar[k] = s[k-1] * (-c[k-2] * beta[k]) - c[k-1] * alpha[k]
699                 *           = s[k-1] * dbar[k] - c[k-1] * alpha[k],
700                 *   delta[k] = c[k-1] * dbar[k] + s[k-1] * alpha[k],
701                 * are not stated in Paige and Saunders (1975), but can be retrieved
702                 * by expanding the (k, k-1) and (k, k) coefficients of the matrix in
703                 * equation (5.5).
704                 */
705                final double deltak = c * dbar + s * alpha;
706                gbar = s * dbar - c * alpha;
707                final double eps = s * beta;
708                dbar = -c * beta;
709                final double zeta = gammaZeta / gamma;
710                /*
711                 * At this point
712                 *   gbar   = gbar[k]
713                 *   deltak = delta[k]
714                 *   eps    = eps[k+1]
715                 *   dbar   = dbar[k+1]
716                 *   zeta   = zeta[k-1]
717                 */
718                final double zetaC = zeta * c;
719                final double zetaS = zeta * s;
720                final int n = xL.getDimension();
721                for (int i = 0; i < n; i++) {
722                    final double xi = xL.getEntry(i);
723                    final double vi = v.getEntry(i);
724                    final double wi = wbar.getEntry(i);
725                    xL.setEntry(i, xi + wi * zetaC + vi * zetaS);
726                    wbar.setEntry(i, wi * s - vi * c);
727                }
728                /*
729                 * At this point
730                 *   x = xL[k-1],
731                 *   ptwbar = P' wbar[k],
732                 * see Paige and Saunders (1975), equations (5.9) and (5.10).
733                 */
734                bstep += snprod * c * zeta;
735                snprod *= s;
736                gmax = FastMath.max(gmax, gamma);
737                gmin = FastMath.min(gmin, gamma);
738                ynorm2 += zeta * zeta;
739                gammaZeta = minusEpsZeta - deltak * zeta;
740                minusEpsZeta = -eps * zeta;
741                /*
742                 * At this point
743                 *   snprod       = s[1] * ... * s[k-1],
744                 *   gmax         = max(|alpha[1]|, gamma[1], ..., gamma[k-1]),
745                 *   gmin         = min(|alpha[1]|, gamma[1], ..., gamma[k-1]),
746                 *   ynorm2       = zeta[1]^2 + ... + zeta[k-1]^2,
747                 *   gammaZeta    = gamma[k] * zeta[k],
748                 *   minusEpsZeta = -eps[k+1] * zeta[k-1].
749                 * The relation for gammaZeta can be retrieved from Paige and
750                 * Saunders (1975), equation (5.4a), last line of the vector
751                 * gbar[k] * zbar[k] = -eps[k] * zeta[k-2] - delta[k] * zeta[k-1].
752                 */
753                updateNorms();
754            }
755    
756            /**
757             * Computes the norms of the residuals, and checks for convergence.
758             * Updates {@link #lqnorm} and {@link #cgnorm}.
759             */
760            private void updateNorms() {
761                final double anorm = FastMath.sqrt(tnorm);
762                final double ynorm = FastMath.sqrt(ynorm2);
763                final double epsa = anorm * MACH_PREC;
764                final double epsx = anorm * ynorm * MACH_PREC;
765                final double epsr = anorm * ynorm * delta;
766                final double diag = gbar == 0. ? epsa : gbar;
767                lqnorm = FastMath.sqrt(gammaZeta * gammaZeta +
768                                       minusEpsZeta * minusEpsZeta);
769                final double qrnorm = snprod * beta1;
770                cgnorm = qrnorm * beta / FastMath.abs(diag);
771    
772                /*
773                 * Estimate cond(A). In this version we look at the diagonals of L
774                 * in the factorization of the tridiagonal matrix, T = L * Q.
775                 * Sometimes, T[k] can be misleadingly ill-conditioned when T[k+1]
776                 * is not, so we must be careful not to overestimate acond.
777                 */
778                final double acond;
779                if (lqnorm <= cgnorm) {
780                    acond = gmax / gmin;
781                } else {
782                    acond = gmax / FastMath.min(gmin, FastMath.abs(diag));
783                }
784                if (acond * MACH_PREC >= 0.1) {
785                    throw new IllConditionedOperatorException(acond);
786                }
787                if (beta1 <= epsx) {
788                    /*
789                     * x has converged to an eigenvector of A corresponding to the
790                     * eigenvalue shift.
791                     */
792                    throw new SingularOperatorException();
793                }
794                rnorm = FastMath.min(cgnorm, lqnorm);
795                hasConverged = (cgnorm <= epsx) || (cgnorm <= epsr);
796            }
797    
798            /**
799             * Returns {@code true} if the default stopping criterion is fulfilled.
800             *
801             * @return {@code true} if convergence of the iterations has occured
802             */
803            boolean hasConverged() {
804                return hasConverged;
805            }
806    
807            /**
808             * Returns {@code true} if the right-hand side vector is zero exactly.
809             *
810             * @return the boolean value of {@code b == 0}
811             */
812            boolean bEqualsNullVector() {
813                return bIsNull;
814            }
815    
816            /**
817             * Returns {@code true} if {@code beta} is essentially zero. This method
818             * is used to check for early stop of the iterations.
819             *
820             * @return {@code true} if {@code beta < }{@link #MACH_PREC}
821             */
822            boolean betaEqualsZero() {
823                return beta < MACH_PREC;
824            }
825    
826            /**
827             * Returns the norm of the updated, preconditioned residual.
828             *
829             * @return the norm of the residual, ||P * r||
830             */
831            double getNormOfResidual() {
832                return rnorm;
833            }
834        }
835    
836        /** Key for the exception context. */
837        private static final String OPERATOR = "operator";
838    
839        /** Key for the exception context. */
840        private static final String THRESHOLD = "threshold";
841    
842        /** Key for the exception context. */
843        private static final String VECTOR = "vector";
844    
845        /** Key for the exception context. */
846        private static final String VECTOR1 = "vector1";
847    
848        /** Key for the exception context. */
849        private static final String VECTOR2 = "vector2";
850    
851        /** {@code true} if symmetry of matrix and conditioner must be checked. */
852        private final boolean check;
853    
854        /**
855         * The value of the custom tolerance &delta; for the default stopping
856         * criterion.
857         */
858        private final double delta;
859    
860        /**
861         * Creates a new instance of this class, with <a href="#stopcrit">default
862         * stopping criterion</a>. Note that setting {@code check} to {@code true}
863         * entails an extra matrix-vector product in the initial phase.
864         *
865         * @param maxIterations the maximum number of iterations
866         * @param delta the &delta; parameter for the default stopping criterion
867         * @param check {@code true} if self-adjointedness of both matrix and
868         * preconditioner should be checked
869         */
870        public SymmLQ(final int maxIterations, final double delta,
871                      final boolean check) {
872            super(maxIterations);
873            this.delta = delta;
874            this.check = check;
875        }
876    
877        /**
878         * Creates a new instance of this class, with <a href="#stopcrit">default
879         * stopping criterion</a> and custom iteration manager. Note that setting
880         * {@code check} to {@code true} entails an extra matrix-vector product in
881         * the initial phase.
882         *
883         * @param manager the custom iteration manager
884         * @param delta the &delta; parameter for the default stopping criterion
885         * @param check {@code true} if self-adjointedness of both matrix and
886         * preconditioner should be checked
887         */
888        public SymmLQ(final IterationManager manager, final double delta,
889                      final boolean check) {
890            super(manager);
891            this.delta = delta;
892            this.check = check;
893        }
894    
895        /**
896         * Returns {@code true} if symmetry of the matrix, and symmetry as well as
897         * positive definiteness of the preconditioner should be checked.
898         *
899         * @return {@code true} if the tests are to be performed
900         */
901        public final boolean getCheck() {
902            return check;
903        }
904    
905        /**
906         * {@inheritDoc}
907         *
908         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
909         * {@code true}, and {@code a} or {@code m} is not self-adjoint
910         * @throws NonPositiveDefiniteOperatorException if {@code m} is not
911         * positive definite
912         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
913         */
914        @Override
915        public RealVector solve(final RealLinearOperator a,
916            final RealLinearOperator m, final RealVector b) throws
917            NullArgumentException, NonSquareOperatorException,
918            DimensionMismatchException, MaxCountExceededException,
919            NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException,
920            IllConditionedOperatorException {
921            MathUtils.checkNotNull(a);
922            final RealVector x = new ArrayRealVector(a.getColumnDimension());
923            return solveInPlace(a, m, b, x, false, 0.);
924        }
925    
926        /**
927         * Returns an estimate of the solution to the linear system (A - shift
928         * &middot; I) &middot; x = b.
929         * <p>
930         * If the solution x is expected to contain a large multiple of {@code b}
931         * (as in Rayleigh-quotient iteration), then better precision may be
932         * achieved with {@code goodb} set to {@code true}; this however requires an
933         * extra call to the preconditioner.
934         * </p>
935         * <p>
936         * {@code shift} should be zero if the system A &middot; x = b is to be
937         * solved. Otherwise, it could be an approximation to an eigenvalue of A,
938         * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
939         * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
940         * sufficiently like an eigenvector corresponding to an eigenvalue near
941         * shift, then the computed x may have very large components. When
942         * normalized, x may be closer to an eigenvector than b.
943         * </p>
944         *
945         * @param a the linear operator A of the system
946         * @param m the preconditioner, M (can be {@code null})
947         * @param b the right-hand side vector
948         * @param goodb usually {@code false}, except if {@code x} is expected to
949         * contain a large multiple of {@code b}
950         * @param shift the amount to be subtracted to all diagonal elements of A
951         * @return a reference to {@code x} (shallow copy)
952         * @throws NullArgumentException if one of the parameters is {@code null}
953         * @throws NonSquareOperatorException if {@code a} or {@code m} is not square
954         * @throws DimensionMismatchException if {@code m} or {@code b} have dimensions
955         * inconsistent with {@code a}
956         * @throws MaxCountExceededException at exhaustion of the iteration count,
957         * unless a custom
958         * {@link org.apache.commons.math3.util.Incrementor.MaxCountExceededCallback callback}
959         * has been set at construction of the {@link IterationManager}
960         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
961         * {@code true}, and {@code a} or {@code m} is not self-adjoint
962         * @throws NonPositiveDefiniteOperatorException if {@code m} is not
963         * positive definite
964         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
965         */
966        public RealVector solve(final RealLinearOperator a,
967            final RealLinearOperator m, final RealVector b, final boolean goodb,
968            final double shift) throws NullArgumentException,
969            NonSquareOperatorException, DimensionMismatchException,
970            MaxCountExceededException, NonSelfAdjointOperatorException,
971            NonPositiveDefiniteOperatorException, IllConditionedOperatorException {
972            MathUtils.checkNotNull(a);
973            final RealVector x = new ArrayRealVector(a.getColumnDimension());
974            return solveInPlace(a, m, b, x, goodb, shift);
975        }
976    
977        /**
978         * {@inheritDoc}
979         *
980         * @param x not meaningful in this implementation; should not be considered
981         * as an initial guess (<a href="#initguess">more</a>)
982         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
983         * {@code true}, and {@code a} or {@code m} is not self-adjoint
984         * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive
985         * definite
986         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
987         */
988        @Override
989        public RealVector solve(final RealLinearOperator a,
990            final RealLinearOperator m, final RealVector b, final RealVector x)
991            throws NullArgumentException, NonSquareOperatorException,
992            DimensionMismatchException, NonSelfAdjointOperatorException,
993            NonPositiveDefiniteOperatorException, IllConditionedOperatorException,
994            MaxCountExceededException {
995            MathUtils.checkNotNull(x);
996            return solveInPlace(a, m, b, x.copy(), false, 0.);
997        }
998    
999        /**
1000         * {@inheritDoc}
1001         *
1002         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
1003         * {@code true}, and {@code a} is not self-adjoint
1004         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
1005         */
1006        @Override
1007        public RealVector solve(final RealLinearOperator a, final RealVector b)
1008            throws NullArgumentException, NonSquareOperatorException,
1009            DimensionMismatchException, NonSelfAdjointOperatorException,
1010            IllConditionedOperatorException, MaxCountExceededException {
1011            MathUtils.checkNotNull(a);
1012            final RealVector x = new ArrayRealVector(a.getColumnDimension());
1013            x.set(0.);
1014            return solveInPlace(a, null, b, x, false, 0.);
1015        }
1016    
1017        /**
1018         * Returns the solution to the system (A - shift &middot; I) &middot; x = b.
1019         * <p>
1020         * If the solution x is expected to contain a large multiple of {@code b}
1021         * (as in Rayleigh-quotient iteration), then better precision may be
1022         * achieved with {@code goodb} set to {@code true}.
1023         * </p>
1024         * <p>
1025         * {@code shift} should be zero if the system A &middot; x = b is to be
1026         * solved. Otherwise, it could be an approximation to an eigenvalue of A,
1027         * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
1028         * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
1029         * sufficiently like an eigenvector corresponding to an eigenvalue near
1030         * shift, then the computed x may have very large components. When
1031         * normalized, x may be closer to an eigenvector than b.
1032         * </p>
1033         *
1034         * @param a the linear operator A of the system
1035         * @param b the right-hand side vector
1036         * @param goodb usually {@code false}, except if {@code x} is expected to
1037         * contain a large multiple of {@code b}
1038         * @param shift the amount to be subtracted to all diagonal elements of A
1039         * @return a reference to {@code x}
1040         * @throws NullArgumentException if one of the parameters is {@code null}
1041         * @throws NonSquareOperatorException if {@code a} is not square
1042         * @throws DimensionMismatchException if {@code b} has dimensions
1043         * inconsistent with {@code a}
1044         * @throws MaxCountExceededException at exhaustion of the iteration count,
1045         * unless a custom
1046         * {@link org.apache.commons.math3.util.Incrementor.MaxCountExceededCallback callback}
1047         * has been set at construction of the {@link IterationManager}
1048         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
1049         * {@code true}, and {@code a} is not self-adjoint
1050         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
1051         */
1052        public RealVector solve(final RealLinearOperator a, final RealVector b,
1053            final boolean goodb, final double shift) throws NullArgumentException,
1054            NonSquareOperatorException, DimensionMismatchException,
1055            NonSelfAdjointOperatorException, IllConditionedOperatorException,
1056            MaxCountExceededException {
1057            MathUtils.checkNotNull(a);
1058            final RealVector x = new ArrayRealVector(a.getColumnDimension());
1059            return solveInPlace(a, null, b, x, goodb, shift);
1060        }
1061    
1062        /**
1063         * {@inheritDoc}
1064         *
1065         * @param x not meaningful in this implementation; should not be considered
1066         * as an initial guess (<a href="#initguess">more</a>)
1067         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
1068         * {@code true}, and {@code a} is not self-adjoint
1069         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
1070         */
1071        @Override
1072        public RealVector solve(final RealLinearOperator a, final RealVector b,
1073            final RealVector x) throws NullArgumentException,
1074            NonSquareOperatorException, DimensionMismatchException,
1075            NonSelfAdjointOperatorException, IllConditionedOperatorException,
1076            MaxCountExceededException {
1077            MathUtils.checkNotNull(x);
1078            return solveInPlace(a, null, b, x.copy(), false, 0.);
1079        }
1080    
1081        /**
1082         * {@inheritDoc}
1083         *
1084         * @param x the vector to be updated with the solution; {@code x} should
1085         * not be considered as an initial guess (<a href="#initguess">more</a>)
1086         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
1087         * {@code true}, and {@code a} or {@code m} is not self-adjoint
1088         * @throws NonPositiveDefiniteOperatorException if {@code m} is not
1089         * positive definite
1090         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
1091         */
1092        @Override
1093        public RealVector solveInPlace(final RealLinearOperator a,
1094            final RealLinearOperator m, final RealVector b, final RealVector x)
1095            throws NullArgumentException, NonSquareOperatorException,
1096            DimensionMismatchException, NonSelfAdjointOperatorException,
1097            NonPositiveDefiniteOperatorException, IllConditionedOperatorException,
1098            MaxCountExceededException {
1099            return solveInPlace(a, m, b, x, false, 0.);
1100        }
1101    
1102        /**
1103         * Returns an estimate of the solution to the linear system (A - shift
1104         * &middot; I) &middot; x = b. The solution is computed in-place.
1105         * <p>
1106         * If the solution x is expected to contain a large multiple of {@code b}
1107         * (as in Rayleigh-quotient iteration), then better precision may be
1108         * achieved with {@code goodb} set to {@code true}; this however requires an
1109         * extra call to the preconditioner.
1110         * </p>
1111         * <p>
1112         * {@code shift} should be zero if the system A &middot; x = b is to be
1113         * solved. Otherwise, it could be an approximation to an eigenvalue of A,
1114         * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
1115         * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
1116         * sufficiently like an eigenvector corresponding to an eigenvalue near
1117         * shift, then the computed x may have very large components. When
1118         * normalized, x may be closer to an eigenvector than b.
1119         * </p>
1120         *
1121         * @param a the linear operator A of the system
1122         * @param m the preconditioner, M (can be {@code null})
1123         * @param b the right-hand side vector
1124         * @param x the vector to be updated with the solution; {@code x} should
1125         * not be considered as an initial guess (<a href="#initguess">more</a>)
1126         * @param goodb usually {@code false}, except if {@code x} is expected to
1127         * contain a large multiple of {@code b}
1128         * @param shift the amount to be subtracted to all diagonal elements of A
1129         * @return a reference to {@code x} (shallow copy).
1130         * @throws NullArgumentException if one of the parameters is {@code null}
1131         * @throws NonSquareOperatorException if {@code a} or {@code m} is not square
1132         * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x}
1133         * have dimensions inconsistent with {@code a}.
1134         * @throws MaxCountExceededException at exhaustion of the iteration count,
1135         * unless a custom
1136         * {@link org.apache.commons.math3.util.Incrementor.MaxCountExceededCallback callback}
1137         * has been set at construction of the {@link IterationManager}
1138         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
1139         * {@code true}, and {@code a} or {@code m} is not self-adjoint
1140         * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive
1141         * definite
1142         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
1143         */
1144        public RealVector solveInPlace(final RealLinearOperator a,
1145            final RealLinearOperator m, final RealVector b,
1146            final RealVector x, final boolean goodb, final double shift)
1147            throws NullArgumentException, NonSquareOperatorException,
1148            DimensionMismatchException, NonSelfAdjointOperatorException,
1149            NonPositiveDefiniteOperatorException, IllConditionedOperatorException,
1150            MaxCountExceededException {
1151            checkParameters(a, m, b, x);
1152    
1153            final IterationManager manager = getIterationManager();
1154            /* Initialization counts as an iteration. */
1155            manager.resetIterationCount();
1156            manager.incrementIterationCount();
1157    
1158            final State state;
1159            state = new State(a, m, b, goodb, shift, delta, check);
1160            state.init();
1161            state.refineSolution(x);
1162            IterativeLinearSolverEvent event;
1163            event = new DefaultIterativeLinearSolverEvent(this,
1164                                                          manager.getIterations(),
1165                                                          x,
1166                                                          b,
1167                                                          state.getNormOfResidual());
1168            if (state.bEqualsNullVector()) {
1169                /* If b = 0 exactly, stop with x = 0. */
1170                manager.fireTerminationEvent(event);
1171                return x;
1172            }
1173            /* Cause termination if beta is essentially zero. */
1174            final boolean earlyStop;
1175            earlyStop = state.betaEqualsZero() || state.hasConverged();
1176            manager.fireInitializationEvent(event);
1177            if (!earlyStop) {
1178                do {
1179                    manager.incrementIterationCount();
1180                    event = new DefaultIterativeLinearSolverEvent(this,
1181                                                                  manager.getIterations(),
1182                                                                  x,
1183                                                                  b,
1184                                                                  state.getNormOfResidual());
1185                    manager.fireIterationStartedEvent(event);
1186                    state.update();
1187                    state.refineSolution(x);
1188                    event = new DefaultIterativeLinearSolverEvent(this,
1189                                                                  manager.getIterations(),
1190                                                                  x,
1191                                                                  b,
1192                                                                  state.getNormOfResidual());
1193                    manager.fireIterationPerformedEvent(event);
1194                } while (!state.hasConverged());
1195            }
1196            event = new DefaultIterativeLinearSolverEvent(this,
1197                                                          manager.getIterations(),
1198                                                          x,
1199                                                          b,
1200                                                          state.getNormOfResidual());
1201            manager.fireTerminationEvent(event);
1202            return x;
1203        }
1204    
1205        /**
1206         * {@inheritDoc}
1207         *
1208         * @param x the vector to be updated with the solution; {@code x} should
1209         * not be considered as an initial guess (<a href="#initguess">more</a>)
1210         * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
1211         * {@code true}, and {@code a} is not self-adjoint
1212         * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
1213         */
1214        @Override
1215        public RealVector solveInPlace(final RealLinearOperator a,
1216            final RealVector b, final RealVector x) throws NullArgumentException,
1217            NonSquareOperatorException, DimensionMismatchException,
1218            NonSelfAdjointOperatorException, IllConditionedOperatorException,
1219            MaxCountExceededException {
1220            return solveInPlace(a, null, b, x, false, 0.);
1221        }
1222    }