/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.apache.commons.math4.legacy.linear;
  
  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  import org.apache.commons.math4.legacy.exception.NullArgumentException;
  import org.apache.commons.math4.legacy.exception.util.ExceptionContext;
  import org.apache.commons.math4.core.jdkmath.JdkMath;
  
  /**
   * <p>
   * Implementation of the SYMMLQ iterative linear solver proposed by <a
   * href="#PAIG1975">Paige and Saunders (1975)</a>. This implementation is
   * largely based on the FORTRAN code by Pr. Michael A. Saunders, available <a
   * href="http://www.stanford.edu/group/SOL/software/symmlq/f77/">here</a>.
   * </p>
   * <p>
   * SYMMLQ is designed to solve the system of linear equations A &middot; x = b
   * where A is an n &times; n self-adjoint linear operator (defined as a
   * {@link RealLinearOperator}), and b is a given vector. The operator A is not
   * required to be positive definite. If A is known to be definite, the method of
   * conjugate gradients might be preferred, since it will require about the same
   * number of iterations as SYMMLQ but slightly less work per iteration.
   * </p>
   * <p>
   * SYMMLQ is designed to solve the system (A - shift &middot; I) &middot; x = b,
   * where shift is a specified scalar value. If shift and b are suitably chosen,
   * the computed vector x may approximate an (unnormalized) eigenvector of A, as
   * in the methods of inverse iteration and/or Rayleigh-quotient iteration.
   * Again, the linear operator (A - shift &middot; I) need not be positive
   * definite (but <em>must</em> be self-adjoint). The work per iteration is very
   * slightly less if shift = 0.
   * </p>
   * <p><b>Preconditioning</b></p>
   * <p>
   * Preconditioning may reduce the number of iterations required. The solver may
   * be provided with a positive definite preconditioner
   * M = P<sup>T</sup> &middot; P
   * that is known to approximate
   * (A - shift &middot; I)<sup>-1</sup> in some sense, where matrix-vector
   * products of the form M &middot; y = x can be computed efficiently. Then
   * SYMMLQ will implicitly solve the system of equations
   * P &middot; (A - shift &middot; I) &middot; P<sup>T</sup> &middot;
   * x<sub>hat</sub> = P &middot; b, i.e.
   * A<sub>hat</sub> &middot; x<sub>hat</sub> = b<sub>hat</sub>,
   * where
   * A<sub>hat</sub> = P &middot; (A - shift &middot; I) &middot; P<sup>T</sup>,
   * b<sub>hat</sub> = P &middot; b,
   * and return the solution
   * x = P<sup>T</sup> &middot; x<sub>hat</sub>.
   * The associated residual is
   * r<sub>hat</sub> = b<sub>hat</sub> - A<sub>hat</sub> &middot; x<sub>hat</sub>
   *                 = P &middot; [b - (A - shift &middot; I) &middot; x]
   *                 = P &middot; r.
   * </p>
   * <p>
   * In the case of preconditioning, the {@link IterativeLinearSolverEvent}s that
   * this solver fires are such that
   * {@link IterativeLinearSolverEvent#getNormOfResidual()} returns the norm of
   * the <em>preconditioned</em>, updated residual, ||P &middot; r||, not the norm
   * of the <em>true</em> residual ||r||.
   * </p>
   * <p><b><a id="stopcrit">Default stopping criterion</a></b></p>
   * <p>
   * A default stopping criterion is implemented. The iterations stop when || rhat
   * || &le; &delta; || Ahat || || xhat ||, where xhat is the current estimate of
   * the solution of the transformed system, rhat the current estimate of the
   * corresponding residual, and &delta; a user-specified tolerance.
   * </p>
   * <p><b>Iteration count</b></p>
   * <p>
   * In the present context, an iteration should be understood as one evaluation
   * of the matrix-vector product A &middot; x. The initialization phase therefore
   * counts as one iteration. If the user requires checks on the symmetry of A,
   * this entails one further matrix-vector product in the initial phase. This
   * further product is <em>not</em> accounted for in the iteration count. In
   * other words, the number of iterations required to reach convergence will be
   * identical, whether checks have been required or not.
   * </p>
   * <p>
   * The present definition of the iteration count differs from that adopted in
   * the original FOTRAN code, where the initialization phase was <em>not</em>
   * taken into account.
   * </p>
  * <p><b><a id="initguess">Initial guess of the solution</a></b></p>
  * <p>
  * The {@code x} parameter in
  * <ul>
  * <li>{@link #solve(RealLinearOperator, RealVector, RealVector)},</li>
  * <li>{@link #solve(RealLinearOperator, RealLinearOperator, RealVector, RealVector)}},</li>
  * <li>{@link #solveInPlace(RealLinearOperator, RealVector, RealVector)},</li>
  * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector)},</li>
  * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector, boolean, double)},</li>
  * </ul>
  * should not be considered as an initial guess, as it is set to zero in the
  * initial phase. If x<sub>0</sub> is known to be a good approximation to x, one
  * should compute r<sub>0</sub> = b - A &middot; x, solve A &middot; dx = r0,
  * and set x = x<sub>0</sub> + dx.
  *
  * <p><b><a id="context">Exception context</a></b></p>
  * <p>
  * Besides standard {@link DimensionMismatchException}, this class might throw
  * {@link NonSelfAdjointOperatorException} if the linear operator or the
  * preconditioner are not symmetric. In this case, the {@link ExceptionContext}
  * provides more information
  * <ul>
  * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
  * <li>key {@code "vector1"} points to the first offending vector, say x,
  * <li>key {@code "vector2"} points to the second offending vector, say y, such
  * that x<sup>T</sup> &middot; L &middot; y &ne; y<sup>T</sup> &middot; L
  * &middot; x (within a certain accuracy).</li>
  * </ul>
  *
  * <p>
  * {@link NonPositiveDefiniteOperatorException} might also be thrown in case the
  * preconditioner is not positive definite. The relevant keys to the
  * {@link ExceptionContext} are
  * <ul>
  * <li>key {@code "operator"}, which points to the offending linear operator,
  * say L,</li>
  * <li>key {@code "vector"}, which points to the offending vector, say x, such
  * that x<sup>T</sup> &middot; L &middot; x &lt; 0.</li>
  * </ul>
  *
  * <p><b>References</b></p>
  * <dl>
  * <dt><a id="PAIG1975">Paige and Saunders (1975)</a></dt>
  * <dd>C. C. Paige and M. A. Saunders, <a
  * href="http://www.stanford.edu/group/SOL/software/symmlq/PS75.pdf"><em>
  * Solution of Sparse Indefinite Systems of Linear Equations</em></a>, SIAM
  * Journal on Numerical Analysis 12(4): 617-629, 1975</dd>
  * </dl>
  *
  * @since 3.0
  */
 public class SymmLQ
     extends PreconditionedIterativeLinearSolver {
 
     /*
      * IMPLEMENTATION NOTES
      * --------------------
      * The implementation follows as closely as possible the notations of Paige
      * and Saunders (1975). Attention must be paid to the fact that some
      * quantities which are relevant to iteration k can only be computed in
      * iteration (k+1). Therefore, minute attention must be paid to the index of
      * each state variable of this algorithm.
      *
      * 1. Preconditioning
      *    ---------------
      * The Lanczos iterations associated with Ahat and bhat read
      *   beta[1] = ||P * b||
      *   v[1] = P * b / beta[1]
      *   beta[k+1] * v[k+1] = Ahat * v[k] - alpha[k] * v[k] - beta[k] * v[k-1]
      *                      = P * (A - shift * I) * P' * v[k] - alpha[k] * v[k]
      *                        - beta[k] * v[k-1]
      * Multiplying both sides by P', we get
      *   beta[k+1] * (P' * v)[k+1] = M * (A - shift * I) * (P' * v)[k]
      *                               - alpha[k] * (P' * v)[k]
      *                               - beta[k] * (P' * v[k-1]),
      * and
      *   alpha[k+1] = v[k+1]' * Ahat * v[k+1]
      *              = v[k+1]' * P * (A - shift * I) * P' * v[k+1]
      *              = (P' * v)[k+1]' * (A - shift * I) * (P' * v)[k+1].
      *
      * In other words, the Lanczos iterations are unchanged, except for the fact
      * that we really compute (P' * v) instead of v. It can easily be checked
      * that all other formulas are unchanged. It must be noted that P is never
      * explicitly used, only matrix-vector products involving are invoked.
      *
      * 2. Accounting for the shift parameter
      *    ----------------------------------
      * Is trivial: each time A.operate(x) is invoked, one must subtract shift * x
      * to the result.
      *
      * 3. Accounting for the goodb flag
      *    -----------------------------
      * When goodb is set to true, the component of xL along b is computed
      * separately. From Paige and Saunders (1975), equation (5.9), we have
      *   wbar[k+1] = s[k] * wbar[k] - c[k] * v[k+1],
      *   wbar[1] = v[1].
      * Introducing wbar2[k] = wbar[k] - s[1] * ... * s[k-1] * v[1], it can
      * easily be verified by induction that wbar2 follows the same recursive
      * relation
      *   wbar2[k+1] = s[k] * wbar2[k] - c[k] * v[k+1],
      *   wbar2[1] = 0,
      * and we then have
      *   w[k] = c[k] * wbar2[k] + s[k] * v[k+1]
      *          + s[1] * ... * s[k-1] * c[k] * v[1].
      * Introducing w2[k] = w[k] - s[1] * ... * s[k-1] * c[k] * v[1], we find,
      * from (5.10)
      *   xL[k] = zeta[1] * w[1] + ... + zeta[k] * w[k]
      *         = zeta[1] * w2[1] + ... + zeta[k] * w2[k]
      *           + (s[1] * c[2] * zeta[2] + ...
      *           + s[1] * ... * s[k-1] * c[k] * zeta[k]) * v[1]
      *         = xL2[k] + bstep[k] * v[1],
      * where xL2[k] is defined by
      *   xL2[0] = 0,
      *   xL2[k+1] = xL2[k] + zeta[k+1] * w2[k+1],
      * and bstep is defined by
      *   bstep[1] = 0,
      *   bstep[k] = bstep[k-1] + s[1] * ... * s[k-1] * c[k] * zeta[k].
      * We also have, from (5.11)
      *   xC[k] = xL[k-1] + zbar[k] * wbar[k]
      *         = xL2[k-1] + zbar[k] * wbar2[k]
      *           + (bstep[k-1] + s[1] * ... * s[k-1] * zbar[k]) * v[1].
      */
 
     /**
      * <p>
      * A simple container holding the non-final variables used in the
      * iterations. Making the current state of the solver visible from the
      * outside is necessary, because during the iterations, {@code x} does not
      * <em>exactly</em> hold the current estimate of the solution. Indeed,
      * {@code x} needs in general to be moved from the LQ point to the CG point.
      * Besides, additional upudates must be carried out in case {@code goodb} is
      * set to {@code true}.
      * </p>
      * <p>
      * In all subsequent comments, the description of the state variables refer
      * to their value after a call to {@link #update()}. In these comments, k is
      * the current number of evaluations of matrix-vector products.
      * </p>
      */
     private static final class State {
         /** The cubic root of {@link #MACH_PREC}. */
         static final double CBRT_MACH_PREC;
 
         /** The machine precision. */
         static final double MACH_PREC;
 
         /** Reference to the linear operator. */
         private final RealLinearOperator a;
 
         /** Reference to the right-hand side vector. */
         private final RealVector b;
 
         /** {@code true} if symmetry of matrix and conditioner must be checked. */
         private final boolean check;
 
         /**
          * The value of the custom tolerance &delta; for the default stopping
          * criterion.
          */
         private final double delta;
 
         /** The value of beta[k+1]. */
         private double beta;
 
         /** The value of beta[1]. */
         private double beta1;
 
         /** The value of bstep[k-1]. */
         private double bstep;
 
         /** The estimate of the norm of P * rC[k]. */
         private double cgnorm;
 
         /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */
         private double dbar;
 
         /**
          * The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the
          * initial code.
          */
         private double gammaZeta;
 
         /** The value of gbar[k]. */
         private double gbar;
 
         /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
         private double gmax;
 
         /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
         private double gmin;
 
         /** Copy of the {@code goodb} parameter. */
         private final boolean goodb;
 
         /** {@code true} if the default convergence criterion is verified. */
         private boolean hasConverged;
 
         /** The estimate of the norm of P * rL[k-1]. */
         private double lqnorm;
 
         /** Reference to the preconditioner, M. */
         private final RealLinearOperator m;
 
         /**
          * The value of (-eps[k+1] * zeta[k-1]). Was called {@code rhs2} in the
          * initial code.
          */
         private double minusEpsZeta;
 
         /** The value of M * b. */
         private final RealVector mb;
 
         /** The value of beta[k]. */
         private double oldb;
 
         /** The value of beta[k] * M^(-1) * P' * v[k]. */
         private RealVector r1;
 
         /** The value of beta[k+1] * M^(-1) * P' * v[k+1]. */
         private RealVector r2;
 
         /**
          * The value of the updated, preconditioned residual P * r. This value is
          * given by {@code min(}{@link #cgnorm}{@code , }{@link #lqnorm}{@code )}.
          */
         private double rnorm;
 
         /** Copy of the {@code shift} parameter. */
         private final double shift;
 
         /** The value of s[1] * ... * s[k-1]. */
         private double snprod;
 
         /**
          * An estimate of the square of the norm of A * V[k], based on Paige and
          * Saunders (1975), equation (3.3).
          */
         private double tnorm;
 
         /**
          * The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] *
          * v[1]) if {@code goodb} is {@code true}. Was called {@code w} in the
          * initial code.
          */
         private RealVector wbar;
 
         /**
          * A reference to the vector to be updated with the solution. Contains
          * the value of xL[k-1] if {@code goodb} is {@code false}, (xL[k-1] -
          * bstep[k-1] * v[1]) otherwise.
          */
         private final RealVector xL;
 
         /** The value of beta[k+1] * P' * v[k+1]. */
         private RealVector y;
 
         /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */
         private double ynorm2;
 
         /** The value of {@code b == 0} (exact floating-point equality). */
         private boolean bIsNull;
 
         static {
             MACH_PREC = JdkMath.ulp(1.);
             CBRT_MACH_PREC = JdkMath.cbrt(MACH_PREC);
         }
 
         /**
          * Creates and inits to k = 1 a new instance of this class.
          *
          * @param a the linear operator A of the system
          * @param m the preconditioner, M (can be {@code null})
          * @param b the right-hand side vector
          * @param goodb usually {@code false}, except if {@code x} is expected
          * to contain a large multiple of {@code b}
          * @param shift the amount to be subtracted to all diagonal elements of
          * A
          * @param delta the &delta; parameter for the default stopping criterion
          * @param check {@code true} if self-adjointedness of both matrix and
          * preconditioner should be checked
          */
         State(final RealLinearOperator a,
             final RealLinearOperator m,
             final RealVector b,
             final boolean goodb,
             final double shift,
             final double delta,
             final boolean check) {
             this.a = a;
             this.m = m;
             this.b = b;
             this.xL = new ArrayRealVector(b.getDimension());
             this.goodb = goodb;
             this.shift = shift;
             this.mb = m == null ? b : m.operate(b);
             this.hasConverged = false;
             this.check = check;
             this.delta = delta;
         }
 
         /**
          * Performs a symmetry check on the specified linear operator, and throws an
          * exception in case this check fails. Given a linear operator L, and a
          * vector x, this method checks that
          * x' &middot; L &middot; y = y' &middot; L &middot; x
          * (within a given accuracy), where y = L &middot; x.
          *
          * @param l the linear operator L
          * @param x the candidate vector x
          * @param y the candidate vector y = L &middot; x
          * @param z the vector z = L &middot; y
          * @throws NonSelfAdjointOperatorException when the test fails
          */
         private static void checkSymmetry(final RealLinearOperator l,
             final RealVector x, final RealVector y, final RealVector z)
             throws NonSelfAdjointOperatorException {
             final double s = y.dotProduct(y);
             final double t = x.dotProduct(z);
             final double epsa = (s + MACH_PREC) * CBRT_MACH_PREC;
             if (JdkMath.abs(s - t) > epsa) {
                 final NonSelfAdjointOperatorException e;
                 e = new NonSelfAdjointOperatorException();
                 final ExceptionContext context = e.getContext();
                 context.setValue(SymmLQ.OPERATOR, l);
                 context.setValue(SymmLQ.VECTOR1, x);
                 context.setValue(SymmLQ.VECTOR2, y);
                 context.setValue(SymmLQ.THRESHOLD, Double.valueOf(epsa));
                 throw e;
             }
         }
 
         /**
          * Throws a new {@link NonPositiveDefiniteOperatorException} with
          * appropriate context.
          *
          * @param l the offending linear operator
          * @param v the offending vector
          * @throws NonPositiveDefiniteOperatorException in any circumstances
          */
         private static void throwNPDLOException(final RealLinearOperator l,
             final RealVector v) throws NonPositiveDefiniteOperatorException {
             final NonPositiveDefiniteOperatorException e;
             e = new NonPositiveDefiniteOperatorException();
             final ExceptionContext context = e.getContext();
             context.setValue(OPERATOR, l);
             context.setValue(VECTOR, v);
             throw e;
         }
 
         /**
          * A clone of the BLAS {@code DAXPY} function, which carries out the
          * operation y &larr; a &middot; x + y. This is for internal use only: no
          * dimension checks are provided.
          *
          * @param a the scalar by which {@code x} is to be multiplied
          * @param x the vector to be added to {@code y}
          * @param y the vector to be incremented
          */
         private static void daxpy(final double a, final RealVector x,
             final RealVector y) {
             final int n = x.getDimension();
             for (int i = 0; i < n; i++) {
                 y.setEntry(i, a * x.getEntry(i) + y.getEntry(i));
             }
         }
 
         /**
          * A BLAS-like function, for the operation z &larr; a &middot; x + b
          * &middot; y + z. This is for internal use only: no dimension checks are
          * provided.
          *
          * @param a the scalar by which {@code x} is to be multiplied
          * @param x the first vector to be added to {@code z}
          * @param b the scalar by which {@code y} is to be multiplied
          * @param y the second vector to be added to {@code z}
          * @param z the vector to be incremented
          */
         private static void daxpbypz(final double a, final RealVector x,
             final double b, final RealVector y, final RealVector z) {
             final int n = z.getDimension();
             for (int i = 0; i < n; i++) {
                 final double zi;
                 zi = a * x.getEntry(i) + b * y.getEntry(i) + z.getEntry(i);
                 z.setEntry(i, zi);
             }
         }
 
         /**
          * <p>
          * Move to the CG point if it seems better. In this version of SYMMLQ,
          * the convergence tests involve only cgnorm, so we're unlikely to stop
          * at an LQ point, except if the iteration limit interferes.
          * </p>
          * <p>
          * Additional upudates are also carried out in case {@code goodb} is set
          * to {@code true}.
          * </p>
          *
          * @param x the vector to be updated with the refined value of xL
          */
          void refineSolution(final RealVector x) {
             final int n = this.xL.getDimension();
             if (lqnorm < cgnorm) {
                 if (!goodb) {
                     x.setSubVector(0, this.xL);
                 } else {
                     final double step = bstep / beta1;
                     for (int i = 0; i < n; i++) {
                         final double bi = mb.getEntry(i);
                         final double xi = this.xL.getEntry(i);
                         x.setEntry(i, xi + step * bi);
                     }
                 }
             } else {
                 final double anorm = JdkMath.sqrt(tnorm);
                 final double diag = gbar == 0. ? anorm * MACH_PREC : gbar;
                 final double zbar = gammaZeta / diag;
                 final double step = (bstep + snprod * zbar) / beta1;
                 // ynorm = JdkMath.sqrt(ynorm2 + zbar * zbar);
                 if (!goodb) {
                     for (int i = 0; i < n; i++) {
                         final double xi = this.xL.getEntry(i);
                         final double wi = wbar.getEntry(i);
                         x.setEntry(i, xi + zbar * wi);
                     }
                 } else {
                     for (int i = 0; i < n; i++) {
                         final double xi = this.xL.getEntry(i);
                         final double wi = wbar.getEntry(i);
                         final double bi = mb.getEntry(i);
                         x.setEntry(i, xi + zbar * wi + step * bi);
                     }
                 }
             }
         }
 
         /**
          * Performs the initial phase of the SYMMLQ algorithm. On return, the
          * value of the state variables of {@code this} object correspond to k =
          * 1.
          */
          void init() {
             this.xL.set(0.);
             /*
              * Set up y for the first Lanczos vector. y and beta1 will be zero
              * if b = 0.
              */
             this.r1 = this.b.copy();
             this.y = this.m == null ? this.b.copy() : this.m.operate(this.r1);
             if (this.m != null && this.check) {
                 checkSymmetry(this.m, this.r1, this.y, this.m.operate(this.y));
             }
 
             this.beta1 = this.r1.dotProduct(this.y);
             if (this.beta1 < 0.) {
                 throwNPDLOException(this.m, this.y);
             }
             if (this.beta1 == 0.) {
                 /* If b = 0 exactly, stop with x = 0. */
                 this.bIsNull = true;
                 return;
             }
             this.bIsNull = false;
             this.beta1 = JdkMath.sqrt(this.beta1);
             /* At this point
              *   r1 = b,
              *   y = M * b,
              *   beta1 = beta[1].
              */
             final RealVector v = this.y.mapMultiply(1. / this.beta1);
             this.y = this.a.operate(v);
             if (this.check) {
                 checkSymmetry(this.a, v, this.y, this.a.operate(this.y));
             }
             /*
              * Set up y for the second Lanczos vector. y and beta will be zero
              * or very small if b is an eigenvector.
              */
             daxpy(-this.shift, v, this.y);
             final double alpha = v.dotProduct(this.y);
             daxpy(-alpha / this.beta1, this.r1, this.y);
             /*
              * At this point
              *   alpha = alpha[1]
              *   y     = beta[2] * M^(-1) * P' * v[2]
              */
             /* Make sure r2 will be orthogonal to the first v. */
             final double vty = v.dotProduct(this.y);
             final double vtv = v.dotProduct(v);
             daxpy(-vty / vtv, v, this.y);
             this.r2 = this.y.copy();
             if (this.m != null) {
                 this.y = this.m.operate(this.r2);
             }
             this.oldb = this.beta1;
             this.beta = this.r2.dotProduct(this.y);
             if (this.beta < 0.) {
                 throwNPDLOException(this.m, this.y);
             }
             this.beta = JdkMath.sqrt(this.beta);
             /*
              * At this point
              *   oldb = beta[1]
              *   beta = beta[2]
              *   y  = beta[2] * P' * v[2]
              *   r2 = beta[2] * M^(-1) * P' * v[2]
              */
             this.cgnorm = this.beta1;
             this.gbar = alpha;
             this.dbar = this.beta;
             this.gammaZeta = this.beta1;
             this.minusEpsZeta = 0.;
             this.bstep = 0.;
             this.snprod = 1.;
             this.tnorm = alpha * alpha + this.beta * this.beta;
             this.ynorm2 = 0.;
             this.gmax = JdkMath.abs(alpha) + MACH_PREC;
             this.gmin = this.gmax;
 
             if (this.goodb) {
                 this.wbar = new ArrayRealVector(this.a.getRowDimension());
                 this.wbar.set(0.);
             } else {
                 this.wbar = v;
             }
             updateNorms();
         }
 
         /**
          * Performs the next iteration of the algorithm. The iteration count
          * should be incremented prior to calling this method. On return, the
          * value of the state variables of {@code this} object correspond to the
          * current iteration count {@code k}.
          */
         void update() {
             final RealVector v = y.mapMultiply(1. / beta);
             y = a.operate(v);
             daxpbypz(-shift, v, -beta / oldb, r1, y);
             final double alpha = v.dotProduct(y);
             /*
              * At this point
              *   v     = P' * v[k],
              *   y     = (A - shift * I) * P' * v[k] - beta[k] * M^(-1) * P' * v[k-1],
              *   alpha = v'[k] * P * (A - shift * I) * P' * v[k]
              *           - beta[k] * v[k]' * P * M^(-1) * P' * v[k-1]
              *         = v'[k] * P * (A - shift * I) * P' * v[k]
              *           - beta[k] * v[k]' * v[k-1]
              *         = alpha[k].
              */
             daxpy(-alpha / beta, r2, y);
             /*
              * At this point
              *   y = (A - shift * I) * P' * v[k] - alpha[k] * M^(-1) * P' * v[k]
              *       - beta[k] * M^(-1) * P' * v[k-1]
              *     = M^(-1) * P' * (P * (A - shift * I) * P' * v[k] -alpha[k] * v[k]
              *       - beta[k] * v[k-1])
              *     = beta[k+1] * M^(-1) * P' * v[k+1],
              * from Paige and Saunders (1975), equation (3.2).
              *
              * WATCH-IT: the two following lines work only because y is no longer
              * updated up to the end of the present iteration, and is
              * reinitialized at the beginning of the next iteration.
              */
             r1 = r2;
             r2 = y;
             if (m != null) {
                 y = m.operate(r2);
             }
             oldb = beta;
             beta = r2.dotProduct(y);
             if (beta < 0.) {
                 throwNPDLOException(m, y);
             }
             beta = JdkMath.sqrt(beta);
             /*
              * At this point
              *   r1 = beta[k] * M^(-1) * P' * v[k],
              *   r2 = beta[k+1] * M^(-1) * P' * v[k+1],
              *   y  = beta[k+1] * P' * v[k+1],
              *   oldb = beta[k],
              *   beta = beta[k+1].
              */
             tnorm += alpha * alpha + oldb * oldb + beta * beta;
             /*
              * Compute the next plane rotation for Q. See Paige and Saunders
              * (1975), equation (5.6), with
              *   gamma = gamma[k-1],
              *   c     = c[k-1],
              *   s     = s[k-1].
              */
             final double gamma = JdkMath.sqrt(gbar * gbar + oldb * oldb);
             final double c = gbar / gamma;
             final double s = oldb / gamma;
             /*
              * The relations
              *   gbar[k] = s[k-1] * (-c[k-2] * beta[k]) - c[k-1] * alpha[k]
              *           = s[k-1] * dbar[k] - c[k-1] * alpha[k],
              *   delta[k] = c[k-1] * dbar[k] + s[k-1] * alpha[k],
              * are not stated in Paige and Saunders (1975), but can be retrieved
              * by expanding the (k, k-1) and (k, k) coefficients of the matrix in
              * equation (5.5).
              */
             final double deltak = c * dbar + s * alpha;
             gbar = s * dbar - c * alpha;
             final double eps = s * beta;
             dbar = -c * beta;
             final double zeta = gammaZeta / gamma;
             /*
              * At this point
              *   gbar   = gbar[k]
              *   deltak = delta[k]
              *   eps    = eps[k+1]
              *   dbar   = dbar[k+1]
              *   zeta   = zeta[k-1]
              */
             final double zetaC = zeta * c;
             final double zetaS = zeta * s;
             final int n = xL.getDimension();
             for (int i = 0; i < n; i++) {
                 final double xi = xL.getEntry(i);
                 final double vi = v.getEntry(i);
                 final double wi = wbar.getEntry(i);
                 xL.setEntry(i, xi + wi * zetaC + vi * zetaS);
                 wbar.setEntry(i, wi * s - vi * c);
             }
             /*
              * At this point
              *   x = xL[k-1],
              *   ptwbar = P' wbar[k],
              * see Paige and Saunders (1975), equations (5.9) and (5.10).
              */
             bstep += snprod * c * zeta;
             snprod *= s;
             gmax = JdkMath.max(gmax, gamma);
             gmin = JdkMath.min(gmin, gamma);
             ynorm2 += zeta * zeta;
             gammaZeta = minusEpsZeta - deltak * zeta;
             minusEpsZeta = -eps * zeta;
             /*
              * At this point
              *   snprod       = s[1] * ... * s[k-1],
              *   gmax         = max(|alpha[1]|, gamma[1], ..., gamma[k-1]),
              *   gmin         = min(|alpha[1]|, gamma[1], ..., gamma[k-1]),
              *   ynorm2       = zeta[1]^2 + ... + zeta[k-1]^2,
              *   gammaZeta    = gamma[k] * zeta[k],
              *   minusEpsZeta = -eps[k+1] * zeta[k-1].
              * The relation for gammaZeta can be retrieved from Paige and
              * Saunders (1975), equation (5.4a), last line of the vector
              * gbar[k] * zbar[k] = -eps[k] * zeta[k-2] - delta[k] * zeta[k-1].
              */
             updateNorms();
         }
 
         /**
          * Computes the norms of the residuals, and checks for convergence.
          * Updates {@link #lqnorm} and {@link #cgnorm}.
          */
         private void updateNorms() {
             final double anorm = JdkMath.sqrt(tnorm);
             final double ynorm = JdkMath.sqrt(ynorm2);
             final double epsa = anorm * MACH_PREC;
             final double epsx = anorm * ynorm * MACH_PREC;
             final double epsr = anorm * ynorm * delta;
             final double diag = gbar == 0. ? epsa : gbar;
             lqnorm = JdkMath.sqrt(gammaZeta * gammaZeta +
                                    minusEpsZeta * minusEpsZeta);
             final double qrnorm = snprod * beta1;
             cgnorm = qrnorm * beta / JdkMath.abs(diag);
 
             /*
              * Estimate cond(A). In this version we look at the diagonals of L
              * in the factorization of the tridiagonal matrix, T = L * Q.
              * Sometimes, T[k] can be misleadingly ill-conditioned when T[k+1]
              * is not, so we must be careful not to overestimate acond.
              */
             final double acond;
             if (lqnorm <= cgnorm) {
                 acond = gmax / gmin;
             } else {
                 acond = gmax / JdkMath.min(gmin, JdkMath.abs(diag));
             }
             if (acond * MACH_PREC >= 0.1) {
                 throw new IllConditionedOperatorException(acond);
             }
             if (beta1 <= epsx) {
                 /*
                  * x has converged to an eigenvector of A corresponding to the
                  * eigenvalue shift.
                  */
                 throw new SingularOperatorException();
             }
             rnorm = JdkMath.min(cgnorm, lqnorm);
             hasConverged = cgnorm <= epsx || cgnorm <= epsr;
         }
 
         /**
          * Returns {@code true} if the default stopping criterion is fulfilled.
          *
          * @return {@code true} if convergence of the iterations has occurred
          */
         boolean hasConverged() {
             return hasConverged;
         }
 
         /**
          * Returns {@code true} if the right-hand side vector is zero exactly.
          *
          * @return the boolean value of {@code b == 0}
          */
         boolean bEqualsNullVector() {
             return bIsNull;
         }
 
         /**
          * Returns {@code true} if {@code beta} is essentially zero. This method
          * is used to check for early stop of the iterations.
          *
          * @return {@code true} if {@code beta < }{@link #MACH_PREC}
          */
         boolean betaEqualsZero() {
             return beta < MACH_PREC;
         }
 
         /**
          * Returns the norm of the updated, preconditioned residual.
          *
          * @return the norm of the residual, ||P * r||
          */
         double getNormOfResidual() {
             return rnorm;
         }
     }
 
     /** Key for the exception context. */
     private static final String OPERATOR = "operator";
 
     /** Key for the exception context. */
     private static final String THRESHOLD = "threshold";
 
     /** Key for the exception context. */
     private static final String VECTOR = "vector";
 
     /** Key for the exception context. */
     private static final String VECTOR1 = "vector1";
 
     /** Key for the exception context. */
     private static final String VECTOR2 = "vector2";
 
     /** {@code true} if symmetry of matrix and conditioner must be checked. */
     private final boolean check;
 
     /**
      * The value of the custom tolerance &delta; for the default stopping
      * criterion.
      */
     private final double delta;
 
     /**
      * Creates a new instance of this class, with <a href="#stopcrit">default
      * stopping criterion</a>. Note that setting {@code check} to {@code true}
      * entails an extra matrix-vector product in the initial phase.
      *
      * @param maxIterations the maximum number of iterations
      * @param delta the &delta; parameter for the default stopping criterion
      * @param check {@code true} if self-adjointedness of both matrix and
      * preconditioner should be checked
      */
     public SymmLQ(final int maxIterations, final double delta,
                   final boolean check) {
         super(maxIterations);
         this.delta = delta;
         this.check = check;
     }
 
     /**
      * Creates a new instance of this class, with <a href="#stopcrit">default
      * stopping criterion</a> and custom iteration manager. Note that setting
      * {@code check} to {@code true} entails an extra matrix-vector product in
      * the initial phase.
      *
      * @param manager the custom iteration manager
      * @param delta the &delta; parameter for the default stopping criterion
      * @param check {@code true} if self-adjointedness of both matrix and
      * preconditioner should be checked
      */
     public SymmLQ(final IterationManager manager, final double delta,
                   final boolean check) {
         super(manager);
         this.delta = delta;
         this.check = check;
     }
 
     /**
      * Returns {@code true} if symmetry of the matrix, and symmetry as well as
      * positive definiteness of the preconditioner should be checked.
      *
      * @return {@code true} if the tests are to be performed
      */
     public final boolean getCheck() {
         return check;
     }
 
     /**
      * {@inheritDoc}
      *
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} or {@code m} is not self-adjoint
      * @throws NonPositiveDefiniteOperatorException if {@code m} is not
      * positive definite
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     @Override
     public RealVector solve(final RealLinearOperator a,
         final RealLinearOperator m, final RealVector b) throws
         NullArgumentException, NonSquareOperatorException,
         DimensionMismatchException, MaxCountExceededException,
         NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException,
         IllConditionedOperatorException {
         NullArgumentException.check(a);
         final RealVector x = new ArrayRealVector(a.getColumnDimension());
         return solveInPlace(a, m, b, x, false, 0.);
     }
 
     /**
      * Returns an estimate of the solution to the linear system (A - shift
      * &middot; I) &middot; x = b.
      * <p>
      * If the solution x is expected to contain a large multiple of {@code b}
      * (as in Rayleigh-quotient iteration), then better precision may be
      * achieved with {@code goodb} set to {@code true}; this however requires an
      * extra call to the preconditioner.
      * </p>
      * <p>
      * {@code shift} should be zero if the system A &middot; x = b is to be
      * solved. Otherwise, it could be an approximation to an eigenvalue of A,
      * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
      * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
      * sufficiently like an eigenvector corresponding to an eigenvalue near
      * shift, then the computed x may have very large components. When
      * normalized, x may be closer to an eigenvector than b.
      * </p>
      *
      * @param a the linear operator A of the system
      * @param m the preconditioner, M (can be {@code null})
      * @param b the right-hand side vector
      * @param goodb usually {@code false}, except if {@code x} is expected to
      * contain a large multiple of {@code b}
      * @param shift the amount to be subtracted to all diagonal elements of A
      * @return a reference to {@code x} (shallow copy)
      * @throws NullArgumentException if one of the parameters is {@code null}
      * @throws NonSquareOperatorException if {@code a} or {@code m} is not square
      * @throws DimensionMismatchException if {@code m} or {@code b} have dimensions
      * inconsistent with {@code a}
      * @throws MaxCountExceededException at exhaustion of the iteration count,
      * unless a custom
      * {@link org.apache.commons.math4.legacy.core.IntegerSequence.Incrementor.MaxCountExceededCallback callback}
      * has been set at construction of the {@link IterationManager}
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} or {@code m} is not self-adjoint
      * @throws NonPositiveDefiniteOperatorException if {@code m} is not
      * positive definite
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     public RealVector solve(final RealLinearOperator a,
         final RealLinearOperator m, final RealVector b, final boolean goodb,
         final double shift) throws NullArgumentException,
         NonSquareOperatorException, DimensionMismatchException,
         MaxCountExceededException, NonSelfAdjointOperatorException,
         NonPositiveDefiniteOperatorException, IllConditionedOperatorException {
         NullArgumentException.check(a);
         final RealVector x = new ArrayRealVector(a.getColumnDimension());
         return solveInPlace(a, m, b, x, goodb, shift);
     }
 
     /**
      * {@inheritDoc}
      *
      * @param x not meaningful in this implementation; should not be considered
      * as an initial guess (<a href="#initguess">more</a>)
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} or {@code m} is not self-adjoint
      * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive
      * definite
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     @Override
     public RealVector solve(final RealLinearOperator a,
         final RealLinearOperator m, final RealVector b, final RealVector x)
         throws NullArgumentException, NonSquareOperatorException,
         DimensionMismatchException, NonSelfAdjointOperatorException,
         NonPositiveDefiniteOperatorException, IllConditionedOperatorException,
         MaxCountExceededException {
         NullArgumentException.check(x);
         return solveInPlace(a, m, b, x.copy(), false, 0.);
     }
 
     /**
      * {@inheritDoc}
      *
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} is not self-adjoint
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     @Override
     public RealVector solve(final RealLinearOperator a, final RealVector b)
         throws NullArgumentException, NonSquareOperatorException,
         DimensionMismatchException, NonSelfAdjointOperatorException,
         IllConditionedOperatorException, MaxCountExceededException {
         NullArgumentException.check(a);
         final RealVector x = new ArrayRealVector(a.getColumnDimension());
         x.set(0.);
         return solveInPlace(a, null, b, x, false, 0.);
     }
 
     /**
      * Returns the solution to the system (A - shift &middot; I) &middot; x = b.
      * <p>
      * If the solution x is expected to contain a large multiple of {@code b}
      * (as in Rayleigh-quotient iteration), then better precision may be
      * achieved with {@code goodb} set to {@code true}.
      * </p>
      * <p>
      * {@code shift} should be zero if the system A &middot; x = b is to be
      * solved. Otherwise, it could be an approximation to an eigenvalue of A,
      * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
      * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
      * sufficiently like an eigenvector corresponding to an eigenvalue near
      * shift, then the computed x may have very large components. When
      * normalized, x may be closer to an eigenvector than b.
      * </p>
      *
      * @param a the linear operator A of the system
      * @param b the right-hand side vector
      * @param goodb usually {@code false}, except if {@code x} is expected to
      * contain a large multiple of {@code b}
      * @param shift the amount to be subtracted to all diagonal elements of A
      * @return a reference to {@code x}
      * @throws NullArgumentException if one of the parameters is {@code null}
      * @throws NonSquareOperatorException if {@code a} is not square
      * @throws DimensionMismatchException if {@code b} has dimensions
      * inconsistent with {@code a}
      * @throws MaxCountExceededException at exhaustion of the iteration count,
      * unless a custom
      * {@link org.apache.commons.math4.legacy.core.IntegerSequence.Incrementor.MaxCountExceededCallback callback}
      * has been set at construction of the {@link IterationManager}
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} is not self-adjoint
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     public RealVector solve(final RealLinearOperator a, final RealVector b,
         final boolean goodb, final double shift) throws NullArgumentException,
         NonSquareOperatorException, DimensionMismatchException,
         NonSelfAdjointOperatorException, IllConditionedOperatorException,
         MaxCountExceededException {
         NullArgumentException.check(a);
         final RealVector x = new ArrayRealVector(a.getColumnDimension());
         return solveInPlace(a, null, b, x, goodb, shift);
     }
 
     /**
      * {@inheritDoc}
      *
      * @param x not meaningful in this implementation; should not be considered
      * as an initial guess (<a href="#initguess">more</a>)
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} is not self-adjoint
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     @Override
     public RealVector solve(final RealLinearOperator a, final RealVector b,
         final RealVector x) throws NullArgumentException,
         NonSquareOperatorException, DimensionMismatchException,
         NonSelfAdjointOperatorException, IllConditionedOperatorException,
         MaxCountExceededException {
         NullArgumentException.check(x);
         return solveInPlace(a, null, b, x.copy(), false, 0.);
     }
 
     /**
      * {@inheritDoc}
      *
      * @param x the vector to be updated with the solution; {@code x} should
      * not be considered as an initial guess (<a href="#initguess">more</a>)
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} or {@code m} is not self-adjoint
      * @throws NonPositiveDefiniteOperatorException if {@code m} is not
      * positive definite
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     @Override
     public RealVector solveInPlace(final RealLinearOperator a,
         final RealLinearOperator m, final RealVector b, final RealVector x)
         throws NullArgumentException, NonSquareOperatorException,
         DimensionMismatchException, NonSelfAdjointOperatorException,
         NonPositiveDefiniteOperatorException, IllConditionedOperatorException,
         MaxCountExceededException {
         return solveInPlace(a, m, b, x, false, 0.);
     }
 
     /**
      * Returns an estimate of the solution to the linear system (A - shift
      * &middot; I) &middot; x = b. The solution is computed in-place.
      * <p>
      * If the solution x is expected to contain a large multiple of {@code b}
      * (as in Rayleigh-quotient iteration), then better precision may be
      * achieved with {@code goodb} set to {@code true}; this however requires an
      * extra call to the preconditioner.
      * </p>
      * <p>
      * {@code shift} should be zero if the system A &middot; x = b is to be
      * solved. Otherwise, it could be an approximation to an eigenvalue of A,
      * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
      * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
      * sufficiently like an eigenvector corresponding to an eigenvalue near
      * shift, then the computed x may have very large components. When
      * normalized, x may be closer to an eigenvector than b.
      * </p>
      *
      * @param a the linear operator A of the system
      * @param m the preconditioner, M (can be {@code null})
      * @param b the right-hand side vector
      * @param x the vector to be updated with the solution; {@code x} should
      * not be considered as an initial guess (<a href="#initguess">more</a>)
      * @param goodb usually {@code false}, except if {@code x} is expected to
      * contain a large multiple of {@code b}
      * @param shift the amount to be subtracted to all diagonal elements of A
      * @return a reference to {@code x} (shallow copy).
      * @throws NullArgumentException if one of the parameters is {@code null}
      * @throws NonSquareOperatorException if {@code a} or {@code m} is not square
      * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x}
      * have dimensions inconsistent with {@code a}.
      * @throws MaxCountExceededException at exhaustion of the iteration count,
      * unless a custom
      * {@link org.apache.commons.math4.legacy.core.IntegerSequence.Incrementor.MaxCountExceededCallback callback}
      * has been set at construction of the {@link IterationManager}
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} or {@code m} is not self-adjoint
      * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive
      * definite
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     public RealVector solveInPlace(final RealLinearOperator a,
         final RealLinearOperator m, final RealVector b,
         final RealVector x, final boolean goodb, final double shift)
         throws NullArgumentException, NonSquareOperatorException,
         DimensionMismatchException, NonSelfAdjointOperatorException,
         NonPositiveDefiniteOperatorException, IllConditionedOperatorException,
         MaxCountExceededException {
         checkParameters(a, m, b, x);
 
         final IterationManager manager = getIterationManager();
         /* Initialization counts as an iteration. */
         manager.resetIterationCount();
         manager.incrementIterationCount();
 
         final State state;
         state = new State(a, m, b, goodb, shift, delta, check);
         state.init();
         state.refineSolution(x);
         IterativeLinearSolverEvent event;
         event = new DefaultIterativeLinearSolverEvent(this,
                                                       manager.getIterations(),
                                                       x,
                                                       b,
                                                       state.getNormOfResidual());
         if (state.bEqualsNullVector()) {
             /* If b = 0 exactly, stop with x = 0. */
             manager.fireTerminationEvent(event);
             return x;
         }
         /* Cause termination if beta is essentially zero. */
         final boolean earlyStop;
         earlyStop = state.betaEqualsZero() || state.hasConverged();
         manager.fireInitializationEvent(event);
         if (!earlyStop) {
             do {
                 manager.incrementIterationCount();
                 event = new DefaultIterativeLinearSolverEvent(this,
                                                               manager.getIterations(),
                                                               x,
                                                               b,
                                                               state.getNormOfResidual());
                 manager.fireIterationStartedEvent(event);
                 state.update();
                 state.refineSolution(x);
                 event = new DefaultIterativeLinearSolverEvent(this,
                                                               manager.getIterations(),
                                                               x,
                                                               b,
                                                               state.getNormOfResidual());
                 manager.fireIterationPerformedEvent(event);
             } while (!state.hasConverged());
         }
         event = new DefaultIterativeLinearSolverEvent(this,
                                                       manager.getIterations(),
                                                       x,
                                                       b,
                                                       state.getNormOfResidual());
         manager.fireTerminationEvent(event);
         return x;
     }
 
     /**
      * {@inheritDoc}
      *
      * @param x the vector to be updated with the solution; {@code x} should
      * not be considered as an initial guess (<a href="#initguess">more</a>)
      * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is
      * {@code true}, and {@code a} is not self-adjoint
      * @throws IllConditionedOperatorException if {@code a} is ill-conditioned
      */
     @Override
     public RealVector solveInPlace(final RealLinearOperator a,
         final RealVector b, final RealVector x) throws NullArgumentException,
         NonSquareOperatorException, DimensionMismatchException,
         NonSelfAdjointOperatorException, IllConditionedOperatorException,
         MaxCountExceededException {
         return solveInPlace(a, null, b, x, false, 0.);
     }
 }