    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    }