SymmLQ.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.linear;

  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.exception.MathIllegalStateException;
  27. import org.hipparchus.exception.NullArgumentException;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.IterationManager;
  30. import org.hipparchus.util.MathUtils;

  32. /**
  33.  * <p>
  34.  * Implementation of the SYMMLQ iterative linear solver proposed by <a
  35.  * href="#PAIG1975">Paige and Saunders (1975)</a>. This implementation is
  36.  * largely based on the FORTRAN code by Pr. Michael A. Saunders, available <a
  37.  * href="http://www.stanford.edu/group/SOL/software/symmlq/f77/">here</a>.
  38.  * </p>
  39.  * <p>
  40.  * SYMMLQ is designed to solve the system of linear equations A &middot; x = b
  41.  * where A is an n &times; n self-adjoint linear operator (defined as a
  42.  * {@link RealLinearOperator}), and b is a given vector. The operator A is not
  43.  * required to be positive definite. If A is known to be definite, the method of
  44.  * conjugate gradients might be preferred, since it will require about the same
  45.  * number of iterations as SYMMLQ but slightly less work per iteration.
  46.  * </p>
  47.  * <p>
  48.  * SYMMLQ is designed to solve the system (A - shift &middot; I) &middot; x = b,
  49.  * where shift is a specified scalar value. If shift and b are suitably chosen,
  50.  * the computed vector x may approximate an (unnormalized) eigenvector of A, as
  51.  * in the methods of inverse iteration and/or Rayleigh-quotient iteration.
  52.  * Again, the linear operator (A - shift &middot; I) need not be positive
  53.  * definite (but <em>must</em> be self-adjoint). The work per iteration is very
  54.  * slightly less if shift = 0.
  55.  * </p>
  56.  * <p><strong>Preconditioning</strong></p>
  57.  * <p>
  58.  * Preconditioning may reduce the number of iterations required. The solver may
  59.  * be provided with a positive definite preconditioner
  60.  * M = P<sup>T</sup> &middot; P
  61.  * that is known to approximate
  62.  * (A - shift &middot; I)<sup>-1</sup> in some sense, where matrix-vector
  63.  * products of the form M &middot; y = x can be computed efficiently. Then
  64.  * SYMMLQ will implicitly solve the system of equations
  65.  * P &middot; (A - shift &middot; I) &middot; P<sup>T</sup> &middot;
  66.  * x<sub>hat</sub> = P &middot; b, i.e.
  67.  * A<sub>hat</sub> &middot; x<sub>hat</sub> = b<sub>hat</sub>,
  68.  * where
  69.  * A<sub>hat</sub> = P &middot; (A - shift &middot; I) &middot; P<sup>T</sup>,
  70.  * b<sub>hat</sub> = P &middot; b,
  71.  * and return the solution
  72.  * x = P<sup>T</sup> &middot; x<sub>hat</sub>.
  73.  * The associated residual is
  74.  * r<sub>hat</sub> = b<sub>hat</sub> - A<sub>hat</sub> &middot; x<sub>hat</sub>
  75.  *                 = P &middot; [b - (A - shift &middot; I) &middot; x]
  76.  *                 = P &middot; r.
  77.  * </p>
  78.  * <p>
  79.  * In the case of preconditioning, the {@link IterativeLinearSolverEvent}s that
  80.  * this solver fires are such that
  81.  * {@link IterativeLinearSolverEvent#getNormOfResidual()} returns the norm of
  82.  * the <em>preconditioned</em>, updated residual, ||P &middot; r||, not the norm
  83.  * of the <em>true</em> residual ||r||.
  84.  * </p>
  85.  * <p><strong>Default stopping criterion</strong></p>
  86.  * <p>
  87.  * A default stopping criterion is implemented. The iterations stop when || rhat
  88.  * || &le; &delta; || Ahat || || xhat ||, where xhat is the current estimate of
  89.  * the solution of the transformed system, rhat the current estimate of the
  90.  * corresponding residual, and &delta; a user-specified tolerance.
  91.  * </p>
  92.  * <strong>Iteration count</strong>
  93.  * <p>
  94.  * In the present context, an iteration should be understood as one evaluation
  95.  * of the matrix-vector product A &middot; x. The initialization phase therefore
  96.  * counts as one iteration. If the user requires checks on the symmetry of A,
  97.  * this entails one further matrix-vector product in the initial phase. This
  98.  * further product is <em>not</em> accounted for in the iteration count. In
  99.  * other words, the number of iterations required to reach convergence will be
  100.  * identical, whether checks have been required or not.
  101.  * </p>
  102.  * <p>
  103.  * The present definition of the iteration count differs from that adopted in
  104.  * the original FOTRAN code, where the initialization phase was <em>not</em>
  105.  * taken into account.
  106.  * </p>
  107.  * <strong>Initial guess of the solution</strong>
  108.  * <p>
  109.  * The {@code x} parameter in
  110.  * </p>
  111.  * <ul>
  112.  * <li>{@link #solve(RealLinearOperator, RealVector, RealVector)},</li>
  113.  * <li>{@link #solve(RealLinearOperator, RealLinearOperator, RealVector, RealVector)}},</li>
  114.  * <li>{@link #solveInPlace(RealLinearOperator, RealVector, RealVector)},</li>
  115.  * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector)},</li>
  116.  * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector, boolean, double)},</li>
  117.  * </ul>
  118.  * <p>
  119.  * should not be considered as an initial guess, as it is set to zero in the
  120.  * initial phase. If x<sub>0</sub> is known to be a good approximation to x, one
  121.  * should compute r<sub>0</sub> = b - A &middot; x, solve A &middot; dx = r0,
  122.  * and set x = x<sub>0</sub> + dx.
  123.  * </p>
  124.  * <p><strong>Exception context</strong></p>
  125.  * <p>
  126.  * Besides standard {@link MathIllegalArgumentException}, this class might throw
  127.  * {@link MathIllegalArgumentException} if the linear operator or the
  128.  * preconditioner are not symmetric.
  129.  * </p>
  130.  * <ul>
  131.  * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
  132.  * <li>key {@code "vector1"} points to the first offending vector, say x,
  133.  * <li>key {@code "vector2"} points to the second offending vector, say y, such
  134.  * that x<sup>T</sup> &middot; L &middot; y &ne; y<sup>T</sup> &middot; L
  135.  * &middot; x (within a certain accuracy).</li>
  136.  * </ul>
  137.  * <p>
  138.  * {@link MathIllegalArgumentException} might also be thrown in case the
  139.  * preconditioner is not positive definite.
  140.  * </p>
  141.  * <p><strong>References</strong></p>
  142.  * <dl>
  143.  * <dt><a id="PAIG1975">Paige and Saunders (1975)</a></dt>
  144.  * <dd>C. C. Paige and M. A. Saunders, <a
  145.  * href="http://www.stanford.edu/group/SOL/software/symmlq/PS75.pdf"><em>
  146.  * Solution of Sparse Indefinite Systems of Linear Equations</em></a>, SIAM
  147.  * Journal on Numerical Analysis 12(4): 617-629, 1975</dd>
  148.  * </dl>
  149.  *
  150.  */
  151. public class SymmLQ
  152.     extends PreconditionedIterativeLinearSolver {

  154.     /** {@code true} if symmetry of matrix and conditioner must be checked. */
  155.     private final boolean check;

  157.     /**
  158.      * The value of the custom tolerance &delta; for the default stopping
  159.      * criterion.
  160.      */
  161.     private final double delta;

  163.     /*
  164.      * IMPLEMENTATION NOTES
  165.      * --------------------
  166.      * The implementation follows as closely as possible the notations of Paige
  167.      * and Saunders (1975). Attention must be paid to the fact that some
  168.      * quantities which are relevant to iteration k can only be computed in
  169.      * iteration (k+1). Therefore, minute attention must be paid to the index of
  170.      * each state variable of this algorithm.
  171.      *
  172.      * 1. Preconditioning
  173.      *    ---------------
  174.      * The Lanczos iterations associated with Ahat and bhat read
  175.      *   beta[1] = ||P * b||
  176.      *   v[1] = P * b / beta[1]
  177.      *   beta[k+1] * v[k+1] = Ahat * v[k] - alpha[k] * v[k] - beta[k] * v[k-1]
  178.      *                      = P * (A - shift * I) * P' * v[k] - alpha[k] * v[k]
  179.      *                        - beta[k] * v[k-1]
  180.      * Multiplying both sides by P', we get
  181.      *   beta[k+1] * (P' * v)[k+1] = M * (A - shift * I) * (P' * v)[k]
  182.      *                               - alpha[k] * (P' * v)[k]
  183.      *                               - beta[k] * (P' * v[k-1]),
  184.      * and
  185.      *   alpha[k+1] = v[k+1]' * Ahat * v[k+1]
  186.      *              = v[k+1]' * P * (A - shift * I) * P' * v[k+1]
  187.      *              = (P' * v)[k+1]' * (A - shift * I) * (P' * v)[k+1].
  188.      *
  189.      * In other words, the Lanczos iterations are unchanged, except for the fact
  190.      * that we really compute (P' * v) instead of v. It can easily be checked
  191.      * that all other formulas are unchanged. It must be noted that P is never
  192.      * explicitly used, only matrix-vector products involving are invoked.
  193.      *
  194.      * 2. Accounting for the shift parameter
  195.      *    ----------------------------------
  196.      * Is trivial: each time A.operate(x) is invoked, one must subtract shift * x
  197.      * to the result.
  198.      *
  199.      * 3. Accounting for the goodb flag
  200.      *    -----------------------------
  201.      * When goodb is set to true, the component of xL along b is computed
  202.      * separately. From Paige and Saunders (1975), equation (5.9), we have
  203.      *   wbar[k+1] = s[k] * wbar[k] - c[k] * v[k+1],
  204.      *   wbar[1] = v[1].
  205.      * Introducing wbar2[k] = wbar[k] - s[1] * ... * s[k-1] * v[1], it can
  206.      * easily be verified by induction that wbar2 follows the same recursive
  207.      * relation
  208.      *   wbar2[k+1] = s[k] * wbar2[k] - c[k] * v[k+1],
  209.      *   wbar2[1] = 0,
  210.      * and we then have
  211.      *   w[k] = c[k] * wbar2[k] + s[k] * v[k+1]
  212.      *          + s[1] * ... * s[k-1] * c[k] * v[1].
  213.      * Introducing w2[k] = w[k] - s[1] * ... * s[k-1] * c[k] * v[1], we find,
  214.      * from (5.10)
  215.      *   xL[k] = zeta[1] * w[1] + ... + zeta[k] * w[k]
  216.      *         = zeta[1] * w2[1] + ... + zeta[k] * w2[k]
  217.      *           + (s[1] * c[2] * zeta[2] + ...
  218.      *           + s[1] * ... * s[k-1] * c[k] * zeta[k]) * v[1]
  219.      *         = xL2[k] + bstep[k] * v[1],
  220.      * where xL2[k] is defined by
  221.      *   xL2[0] = 0,
  222.      *   xL2[k+1] = xL2[k] + zeta[k+1] * w2[k+1],
  223.      * and bstep is defined by
  224.      *   bstep[1] = 0,
  225.      *   bstep[k] = bstep[k-1] + s[1] * ... * s[k-1] * c[k] * zeta[k].
  226.      * We also have, from (5.11)
  227.      *   xC[k] = xL[k-1] + zbar[k] * wbar[k]
  228.      *         = xL2[k-1] + zbar[k] * wbar2[k]
  229.      *           + (bstep[k-1] + s[1] * ... * s[k-1] * zbar[k]) * v[1].
  230.      */

  232.     /**
  233.      * <p>
  234.      * A simple container holding the non-final variables used in the
  235.      * iterations. Making the current state of the solver visible from the
  236.      * outside is necessary, because during the iterations, {@code x} does not
  237.      * <em>exactly</em> hold the current estimate of the solution. Indeed,
  238.      * {@code x} needs in general to be moved from the LQ point to the CG point.
  239.      * Besides, additional upudates must be carried out in case {@code goodb} is
  240.      * set to {@code true}.
  241.      * </p>
  242.      * <p>
  243.      * In all subsequent comments, the description of the state variables refer
  244.      * to their value after a call to {@link #update()}. In these comments, k is
  245.      * the current number of evaluations of matrix-vector products.
  246.      * </p>
  247.      */
  248.     private static class State {
  249.         /** The cubic root of {@link #MACH_PREC}. */
  250.         static final double CBRT_MACH_PREC;

  252.         /** The machine precision. */
  253.         static final double MACH_PREC;

  255.         /** Reference to the linear operator. */
  256.         private final RealLinearOperator a;

  258.         /** Reference to the right-hand side vector. */
  259.         private final RealVector b;

  261.         /** {@code true} if symmetry of matrix and conditioner must be checked. */
  262.         private final boolean check;

  264.         /**
  265.          * The value of the custom tolerance &delta; for the default stopping
  266.          * criterion.
  267.          */
  268.         private final double delta;

  270.         /** The value of beta[k+1]. */
  271.         private double beta;

  273.         /** The value of beta[1]. */
  274.         private double beta1;

  276.         /** The value of bstep[k-1]. */
  277.         private double bstep;

  279.         /** The estimate of the norm of P * rC[k]. */
  280.         private double cgnorm;

  282.         /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */
  283.         private double dbar;

  285.         /**
  286.          * The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the
  287.          * initial code.
  288.          */
  289.         private double gammaZeta;

  291.         /** The value of gbar[k]. */
  292.         private double gbar;

  294.         /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
  295.         private double gmax;

  297.         /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
  298.         private double gmin;

  300.         /** Copy of the {@code goodb} parameter. */
  301.         private final boolean goodb;

  303.         /** {@code true} if the default convergence criterion is verified. */
  304.         private boolean hasConverged;

  306.         /** The estimate of the norm of P * rL[k-1]. */
  307.         private double lqnorm;

  309.         /** Reference to the preconditioner, M. */
  310.         private final RealLinearOperator m;

  312.         /**
  313.          * The value of (-eps[k+1] * zeta[k-1]). Was called {@code rhs2} in the
  314.          * initial code.
  315.          */
  316.         private double minusEpsZeta;

  318.         /** The value of M * b. */
  319.         private final RealVector mb;

  321.         /** The value of beta[k]. */
  322.         private double oldb;

  324.         /** The value of beta[k] * M^(-1) * P' * v[k]. */
  325.         private RealVector r1;

  327.         /** The value of beta[k+1] * M^(-1) * P' * v[k+1]. */
  328.         private RealVector r2;

  330.         /**
  331.          * The value of the updated, preconditioned residual P * r. This value is
  332.          * given by {@code min(}{@link #cgnorm}{@code , }{@link #lqnorm}{@code )}.
  333.          */
  334.         private double rnorm;

  336.         /** Copy of the {@code shift} parameter. */
  337.         private final double shift;

  339.         /** The value of s[1] * ... * s[k-1]. */
  340.         private double snprod;

  342.         /**
  343.          * An estimate of the square of the norm of A * V[k], based on Paige and
  344.          * Saunders (1975), equation (3.3).
  345.          */
  346.         private double tnorm;

  348.         /**
  349.          * The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] *
  350.          * v[1]) if {@code goodb} is {@code true}. Was called {@code w} in the
  351.          * initial code.
  352.          */
  353.         private RealVector wbar;

  355.         /**
  356.          * A reference to the vector to be updated with the solution. Contains
  357.          * the value of xL[k-1] if {@code goodb} is {@code false}, (xL[k-1] -
  358.          * bstep[k-1] * v[1]) otherwise.
  359.          */
  360.         private final RealVector xL;

  362.         /** The value of beta[k+1] * P' * v[k+1]. */
  363.         private RealVector y;

  365.         /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */
  366.         private double ynorm2;

  368.         /** The value of {@code b == 0} (exact floating-point equality). */
  369.         private boolean bIsNull;

  371.         static {
  372.             MACH_PREC = FastMath.ulp(1.);
  373.             CBRT_MACH_PREC = FastMath.cbrt(MACH_PREC);
  374.         }

  376.         /**
  377.          * Creates and inits to k = 1 a new instance of this class.
  378.          *
  379.          * @param a the linear operator A of the system
  380.          * @param m the preconditioner, M (can be {@code null})
  381.          * @param b the right-hand side vector
  382.          * @param goodb usually {@code false}, except if {@code x} is expected
  383.          * to contain a large multiple of {@code b}
  384.          * @param shift the amount to be subtracted to all diagonal elements of
  385.          * A
  386.          * @param delta the &delta; parameter for the default stopping criterion
  387.          * @param check {@code true} if self-adjointedness of both matrix and
  388.          * preconditioner should be checked
  389.          */
  390.         State(final RealLinearOperator a,
  391.             final RealLinearOperator m,
  392.             final RealVector b,
  393.             final boolean goodb,
  394.             final double shift,
  395.             final double delta,
  396.             final boolean check) {
  397.             this.a = a;
  398.             this.m = m;
  399.             this.b = b;
  400.             this.xL = new ArrayRealVector(b.getDimension());
  401.             this.goodb = goodb;
  402.             this.shift = shift;
  403.             this.mb = m == null ? b : m.operate(b);
  404.             this.hasConverged = false;
  405.             this.check = check;
  406.             this.delta = delta;
  407.         }

  409.         /**
  410.          * Performs a symmetry check on the specified linear operator, and throws an
  411.          * exception in case this check fails. Given a linear operator L, and a
  412.          * vector x, this method checks that
  413.          * x' &middot; L &middot; y = y' &middot; L &middot; x
  414.          * (within a given accuracy), where y = L &middot; x.
  415.          * @param x the candidate vector x
  416.          * @param y the candidate vector y = L &middot; x
  417.          * @param z the vector z = L &middot; y
  418.          *
  419.          * @throws MathIllegalArgumentException when the test fails
  420.          */
  421.         private static void checkSymmetry(final RealVector x, final RealVector y, final RealVector z)
  422.             throws MathIllegalArgumentException {
  423.             final double s = y.dotProduct(y);
  424.             final double t = x.dotProduct(z);
  425.             final double epsa = (s + MACH_PREC) * CBRT_MACH_PREC;
  426.             if (FastMath.abs(s - t) > epsa) {
  427.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SELF_ADJOINT_OPERATOR);
  428.             }
  429.         }

  431.         /**
  432.          * Throws a new {@link MathIllegalArgumentException} with
  433.          * appropriate context.
  434.          * @throws MathIllegalArgumentException in any circumstances
  435.          */
  436.         private static void throwNPDLOException() throws MathIllegalArgumentException {
  437.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_POSITIVE_DEFINITE_OPERATOR);
  438.         }

  440.         /**
  441.          * A clone of the BLAS {@code DAXPY} function, which carries out the
  442.          * operation y &larr; a &middot; x + y. This is for internal use only: no
  443.          * dimension checks are provided.
  444.          *
  445.          * @param a the scalar by which {@code x} is to be multiplied
  446.          * @param x the vector to be added to {@code y}
  447.          * @param y the vector to be incremented
  448.          */
  449.         private static void daxpy(final double a, final RealVector x,
  450.             final RealVector y) {
  451.             final int n = x.getDimension();
  452.             for (int i = 0; i < n; i++) {
  453.                 y.setEntry(i, a * x.getEntry(i) + y.getEntry(i));
  454.             }
  455.         }

  457.         /**
  458.          * A BLAS-like function, for the operation z &larr; a &middot; x + b
  459.          * &middot; y + z. This is for internal use only: no dimension checks are
  460.          * provided.
  461.          *
  462.          * @param a the scalar by which {@code x} is to be multiplied
  463.          * @param x the first vector to be added to {@code z}
  464.          * @param b the scalar by which {@code y} is to be multiplied
  465.          * @param y the second vector to be added to {@code z}
  466.          * @param z the vector to be incremented
  467.          */
  468.         private static void daxpbypz(final double a, final RealVector x,
  469.             final double b, final RealVector y, final RealVector z) {
  470.             final int n = z.getDimension();
  471.             for (int i = 0; i < n; i++) {
  472.                 final double zi;
  473.                 zi = a * x.getEntry(i) + b * y.getEntry(i) + z.getEntry(i);
  474.                 z.setEntry(i, zi);
  475.             }
  476.         }

  478.         /**
  479.          * <p>
  480.          * Move to the CG point if it seems better. In this version of SYMMLQ,
  481.          * the convergence tests involve only cgnorm, so we're unlikely to stop
  482.          * at an LQ point, except if the iteration limit interferes.
  483.          * </p>
  484.          * <p>
  485.          * Additional upudates are also carried out in case {@code goodb} is set
  486.          * to {@code true}.
  487.          * </p>
  488.          *
  489.          * @param x the vector to be updated with the refined value of xL
  490.          */
  491.          void refineSolution(final RealVector x) {
  492.             final int n = this.xL.getDimension();
  493.             if (lqnorm < cgnorm) {
  494.                 if (!goodb) {
  495.                     x.setSubVector(0, this.xL);
  496.                 } else {
  497.                     final double step = bstep / beta1;
  498.                     for (int i = 0; i < n; i++) {
  499.                         final double bi = mb.getEntry(i);
  500.                         final double xi = this.xL.getEntry(i);
  501.                         x.setEntry(i, xi + step * bi);
  502.                     }
  503.                 }
  504.             } else {
  505.                 final double anorm = FastMath.sqrt(tnorm);
  506.                 final double diag = gbar == 0. ? anorm * MACH_PREC : gbar;
  507.                 final double zbar = gammaZeta / diag;
  508.                 final double step = (bstep + snprod * zbar) / beta1;
  509.                 // ynorm = FastMath.sqrt(ynorm2 + zbar * zbar);
  510.                 if (!goodb) {
  511.                     for (int i = 0; i < n; i++) {
  512.                         final double xi = this.xL.getEntry(i);
  513.                         final double wi = wbar.getEntry(i);
  514.                         x.setEntry(i, xi + zbar * wi);
  515.                     }
  516.                 } else {
  517.                     for (int i = 0; i < n; i++) {
  518.                         final double xi = this.xL.getEntry(i);
  519.                         final double wi = wbar.getEntry(i);
  520.                         final double bi = mb.getEntry(i);
  521.                         x.setEntry(i, xi + zbar * wi + step * bi);
  522.                     }
  523.                 }
  524.             }
  525.         }

  527.         /**
  528.          * Performs the initial phase of the SYMMLQ algorithm. On return, the
  529.          * value of the state variables of {@code this} object correspond to k =
  530.          * 1.
  531.          */
  532.          void init() {
  533.             this.xL.set(0.);
  534.             /*
  535.              * Set up y for the first Lanczos vector. y and beta1 will be zero
  536.              * if b = 0.
  537.              */
  538.             this.r1 = this.b.copy();
  539.             this.y = this.m == null ? this.b.copy() : this.m.operate(this.r1);
  540.             if ((this.m != null) && this.check) {
  541.                 checkSymmetry(this.r1, this.y, this.m.operate(this.y));
  542.             }

  544.             this.beta1 = this.r1.dotProduct(this.y);
  545.             if (this.beta1 < 0.) {
  546.                 throwNPDLOException();
  547.             }
  548.             if (this.beta1 == 0.) {
  549.                 /* If b = 0 exactly, stop with x = 0. */
  550.                 this.bIsNull = true;
  551.                 return;
  552.             }
  553.             this.bIsNull = false;
  554.             this.beta1 = FastMath.sqrt(this.beta1);
  555.             /* At this point
  556.              *   r1 = b,
  557.              *   y = M * b,
  558.              *   beta1 = beta[1].
  559.              */
  560.             final RealVector v = this.y.mapMultiply(1. / this.beta1);
  561.             this.y = this.a.operate(v);
  562.             if (this.check) {
  563.                 checkSymmetry(v, this.y, this.a.operate(this.y));
  564.             }
  565.             /*
  566.              * Set up y for the second Lanczos vector. y and beta will be zero
  567.              * or very small if b is an eigenvector.
  568.              */
  569.             daxpy(-this.shift, v, this.y);
  570.             final double alpha = v.dotProduct(this.y);
  571.             daxpy(-alpha / this.beta1, this.r1, this.y);
  572.             /*
  573.              * At this point
  574.              *   alpha = alpha[1]
  575.              *   y     = beta[2] * M^(-1) * P' * v[2]
  576.              */
  577.             /* Make sure r2 will be orthogonal to the first v. */
  578.             final double vty = v.dotProduct(this.y);
  579.             final double vtv = v.dotProduct(v);
  580.             daxpy(-vty / vtv, v, this.y);
  581.             this.r2 = this.y.copy();
  582.             if (this.m != null) {
  583.                 this.y = this.m.operate(this.r2);
  584.             }
  585.             this.oldb = this.beta1;
  586.             this.beta = this.r2.dotProduct(this.y);
  587.             if (this.beta < 0.) {
  588.                 throwNPDLOException();
  589.             }
  590.             this.beta = FastMath.sqrt(this.beta);
  591.             /*
  592.              * At this point
  593.              *   oldb = beta[1]
  594.              *   beta = beta[2]
  595.              *   y  = beta[2] * P' * v[2]
  596.              *   r2 = beta[2] * M^(-1) * P' * v[2]
  597.              */
  598.             this.cgnorm = this.beta1;
  599.             this.gbar = alpha;
  600.             this.dbar = this.beta;
  601.             this.gammaZeta = this.beta1;
  602.             this.minusEpsZeta = 0.;
  603.             this.bstep = 0.;
  604.             this.snprod = 1.;
  605.             this.tnorm = alpha * alpha + this.beta * this.beta;
  606.             this.ynorm2 = 0.;
  607.             this.gmax = FastMath.abs(alpha) + MACH_PREC;
  608.             this.gmin = this.gmax;

  610.             if (this.goodb) {
  611.                 this.wbar = new ArrayRealVector(this.a.getRowDimension());
  612.                 this.wbar.set(0.);
  613.             } else {
  614.                 this.wbar = v;
  615.             }
  616.             updateNorms();
  617.         }

  619.         /**
  620.          * Performs the next iteration of the algorithm. The iteration count
  621.          * should be incremented prior to calling this method. On return, the
  622.          * value of the state variables of {@code this} object correspond to the
  623.          * current iteration count {@code k}.
  624.          */
  625.         void update() {
  626.             final RealVector v = y.mapMultiply(1. / beta);
  627.             y = a.operate(v);
  628.             daxpbypz(-shift, v, -beta / oldb, r1, y);
  629.             final double alpha = v.dotProduct(y);
  630.             /*
  631.              * At this point
  632.              *   v     = P' * v[k],
  633.              *   y     = (A - shift * I) * P' * v[k] - beta[k] * M^(-1) * P' * v[k-1],
  634.              *   alpha = v'[k] * P * (A - shift * I) * P' * v[k]
  635.              *           - beta[k] * v[k]' * P * M^(-1) * P' * v[k-1]
  636.              *         = v'[k] * P * (A - shift * I) * P' * v[k]
  637.              *           - beta[k] * v[k]' * v[k-1]
  638.              *         = alpha[k].
  639.              */
  640.             daxpy(-alpha / beta, r2, y);
  641.             /*
  642.              * At this point
  643.              *   y = (A - shift * I) * P' * v[k] - alpha[k] * M^(-1) * P' * v[k]
  644.              *       - beta[k] * M^(-1) * P' * v[k-1]
  645.              *     = M^(-1) * P' * (P * (A - shift * I) * P' * v[k] -alpha[k] * v[k]
  646.              *       - beta[k] * v[k-1])
  647.              *     = beta[k+1] * M^(-1) * P' * v[k+1],
  648.              * from Paige and Saunders (1975), equation (3.2).
  649.              *
  650.              * WATCH-IT: the two following lines work only because y is no longer
  651.              * updated up to the end of the present iteration, and is
  652.              * reinitialized at the beginning of the next iteration.
  653.              */
  654.             r1 = r2;
  655.             r2 = y;
  656.             if (m != null) {
  657.                 y = m.operate(r2);
  658.             }
  659.             oldb = beta;
  660.             beta = r2.dotProduct(y);
  661.             if (beta < 0.) {
  662.                 throwNPDLOException();
  663.             }
  664.             beta = FastMath.sqrt(beta);
  665.             /*
  666.              * At this point
  667.              *   r1 = beta[k] * M^(-1) * P' * v[k],
  668.              *   r2 = beta[k+1] * M^(-1) * P' * v[k+1],
  669.              *   y  = beta[k+1] * P' * v[k+1],
  670.              *   oldb = beta[k],
  671.              *   beta = beta[k+1].
  672.              */
  673.             tnorm += alpha * alpha + oldb * oldb + beta * beta;
  674.             /*
  675.              * Compute the next plane rotation for Q. See Paige and Saunders
  676.              * (1975), equation (5.6), with
  677.              *   gamma = gamma[k-1],
  678.              *   c     = c[k-1],
  679.              *   s     = s[k-1].
  680.              */
  681.             final double gamma = FastMath.sqrt(gbar * gbar + oldb * oldb);
  682.             final double c = gbar / gamma;
  683.             final double s = oldb / gamma;
  684.             /*
  685.              * The relations
  686.              *   gbar[k] = s[k-1] * (-c[k-2] * beta[k]) - c[k-1] * alpha[k]
  687.              *           = s[k-1] * dbar[k] - c[k-1] * alpha[k],
  688.              *   delta[k] = c[k-1] * dbar[k] + s[k-1] * alpha[k],
  689.              * are not stated in Paige and Saunders (1975), but can be retrieved
  690.              * by expanding the (k, k-1) and (k, k) coefficients of the matrix in
  691.              * equation (5.5).
  692.              */
  693.             final double deltak = c * dbar + s * alpha;
  694.             gbar = s * dbar - c * alpha;
  695.             final double eps = s * beta;
  696.             dbar = -c * beta;
  697.             final double zeta = gammaZeta / gamma;
  698.             /*
  699.              * At this point
  700.              *   gbar   = gbar[k]
  701.              *   deltak = delta[k]
  702.              *   eps    = eps[k+1]
  703.              *   dbar   = dbar[k+1]
  704.              *   zeta   = zeta[k-1]
  705.              */
  706.             final double zetaC = zeta * c;
  707.             final double zetaS = zeta * s;
  708.             final int n = xL.getDimension();
  709.             for (int i = 0; i < n; i++) {
  710.                 final double xi = xL.getEntry(i);
  711.                 final double vi = v.getEntry(i);
  712.                 final double wi = wbar.getEntry(i);
  713.                 xL.setEntry(i, xi + wi * zetaC + vi * zetaS);
  714.                 wbar.setEntry(i, wi * s - vi * c);
  715.             }
  716.             /*
  717.              * At this point
  718.              *   x = xL[k-1],
  719.              *   ptwbar = P' wbar[k],
  720.              * see Paige and Saunders (1975), equations (5.9) and (5.10).
  721.              */
  722.             bstep += snprod * c * zeta;
  723.             snprod *= s;
  724.             gmax = FastMath.max(gmax, gamma);
  725.             gmin = FastMath.min(gmin, gamma);
  726.             ynorm2 += zeta * zeta;
  727.             gammaZeta = minusEpsZeta - deltak * zeta;
  728.             minusEpsZeta = -eps * zeta;
  729.             /*
  730.              * At this point
  731.              *   snprod       = s[1] * ... * s[k-1],
  732.              *   gmax         = max(|alpha[1]|, gamma[1], ..., gamma[k-1]),
  733.              *   gmin         = min(|alpha[1]|, gamma[1], ..., gamma[k-1]),
  734.              *   ynorm2       = zeta[1]^2 + ... + zeta[k-1]^2,
  735.              *   gammaZeta    = gamma[k] * zeta[k],
  736.              *   minusEpsZeta = -eps[k+1] * zeta[k-1].
  737.              * The relation for gammaZeta can be retrieved from Paige and
  738.              * Saunders (1975), equation (5.4a), last line of the vector
  739.              * gbar[k] * zbar[k] = -eps[k] * zeta[k-2] - delta[k] * zeta[k-1].
  740.              */
  741.             updateNorms();
  742.         }

  744.         /**
  745.          * Computes the norms of the residuals, and checks for convergence.
  746.          * Updates {@link #lqnorm} and {@link #cgnorm}.
  747.          */
  748.         private void updateNorms() {
  749.             final double anorm = FastMath.sqrt(tnorm);
  750.             final double ynorm = FastMath.sqrt(ynorm2);
  751.             final double epsa = anorm * MACH_PREC;
  752.             final double epsx = anorm * ynorm * MACH_PREC;
  753.             final double epsr = anorm * ynorm * delta;
  754.             final double diag = gbar == 0. ? epsa : gbar;
  755.             lqnorm = FastMath.sqrt(gammaZeta * gammaZeta +
  756.                                    minusEpsZeta * minusEpsZeta);
  757.             final double qrnorm = snprod * beta1;
  758.             cgnorm = qrnorm * beta / FastMath.abs(diag);

  760.             /*
  761.              * Estimate cond(A). In this version we look at the diagonals of L
  762.              * in the factorization of the tridiagonal matrix, T = L * Q.
  763.              * Sometimes, T[k] can be misleadingly ill-conditioned when T[k+1]
  764.              * is not, so we must be careful not to overestimate acond.
  765.              */
  766.             final double acond;
  767.             if (lqnorm <= cgnorm) {
  768.                 acond = gmax / gmin;
  769.             } else {
  770.                 acond = gmax / FastMath.min(gmin, FastMath.abs(diag));
  771.             }
  772.             if (acond * MACH_PREC >= 0.1) {
  773.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.ILL_CONDITIONED_OPERATOR, acond);
  774.             }
  775.             if (beta1 <= epsx) {
  776.                 /*
  777.                  * x has converged to an eigenvector of A corresponding to the
  778.                  * eigenvalue shift.
  779.                  */
  780.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_OPERATOR);
  781.             }
  782.             rnorm = FastMath.min(cgnorm, lqnorm);
  783.             hasConverged = (cgnorm <= epsx) || (cgnorm <= epsr);
  784.         }

  786.         /**
  787.          * Returns {@code true} if the default stopping criterion is fulfilled.
  788.          *
  789.          * @return {@code true} if convergence of the iterations has occurred
  790.          */
  791.         boolean hasConverged() {
  792.             return hasConverged;
  793.         }

  795.         /**
  796.          * Returns {@code true} if the right-hand side vector is zero exactly.
  797.          *
  798.          * @return the boolean value of {@code b == 0}
  799.          */
  800.         boolean bEqualsNullVector() {
  801.             return bIsNull;
  802.         }

  804.         /**
  805.          * Returns {@code true} if {@code beta} is essentially zero. This method
  806.          * is used to check for early stop of the iterations.
  807.          *
  808.          * @return {@code true} if {@code beta < }{@link #MACH_PREC}
  809.          */
  810.         boolean betaEqualsZero() {
  811.             return beta < MACH_PREC;
  812.         }

  814.         /**
  815.          * Returns the norm of the updated, preconditioned residual.
  816.          *
  817.          * @return the norm of the residual, ||P * r||
  818.          */
  819.         double getNormOfResidual() {
  820.             return rnorm;
  821.         }
  822.     }

  824.     /**
  825.      * Creates a new instance of this class, with <a href="#stopcrit">default
  826.      * stopping criterion</a>. Note that setting {@code check} to {@code true}
  827.      * entails an extra matrix-vector product in the initial phase.
  828.      *
  829.      * @param maxIterations the maximum number of iterations
  830.      * @param delta the &delta; parameter for the default stopping criterion
  831.      * @param check {@code true} if self-adjointedness of both matrix and
  832.      * preconditioner should be checked
  833.      */
  834.     public SymmLQ(final int maxIterations, final double delta,
  835.                   final boolean check) {
  836.         super(maxIterations);
  837.         this.delta = delta;
  838.         this.check = check;
  839.     }

  841.     /**
  842.      * Creates a new instance of this class, with <a href="#stopcrit">default
  843.      * stopping criterion</a> and custom iteration manager. Note that setting
  844.      * {@code check} to {@code true} entails an extra matrix-vector product in
  845.      * the initial phase.
  846.      *
  847.      * @param manager the custom iteration manager
  848.      * @param delta the &delta; parameter for the default stopping criterion
  849.      * @param check {@code true} if self-adjointedness of both matrix and
  850.      * preconditioner should be checked
  851.      */
  852.     public SymmLQ(final IterationManager manager, final double delta,
  853.                   final boolean check) {
  854.         super(manager);
  855.         this.delta = delta;
  856.         this.check = check;
  857.     }

  859.     /**
  860.      * Returns {@code true} if symmetry of the matrix, and symmetry as well as
  861.      * positive definiteness of the preconditioner should be checked.
  862.      *
  863.      * @return {@code true} if the tests are to be performed
  864.      * @since 1.4
  865.      */
  866.     public final boolean shouldCheck() {
  867.         return check;
  868.     }

  870.     /**
  871.      * {@inheritDoc}
  872.      *
  873.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  874.      * {@code true}, and {@code a} or {@code m} is not self-adjoint
  875.      * @throws MathIllegalArgumentException if {@code m} is not
  876.      * positive definite
  877.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  878.      */
  879.     @Override
  880.     public RealVector solve(final RealLinearOperator a,
  881.         final RealLinearOperator m, final RealVector b) throws
  882.             NullArgumentException, MathIllegalStateException, MathIllegalArgumentException {
  883.         MathUtils.checkNotNull(a);
  884.         final RealVector x = new ArrayRealVector(a.getColumnDimension());
  885.         return solveInPlace(a, m, b, x, false, 0.);
  886.     }

  888.     /**
  889.      * Returns an estimate of the solution to the linear system (A - shift
  890.      * &middot; I) &middot; x = b.
  891.      * <p>
  892.      * If the solution x is expected to contain a large multiple of {@code b}
  893.      * (as in Rayleigh-quotient iteration), then better precision may be
  894.      * achieved with {@code goodb} set to {@code true}; this however requires an
  895.      * extra call to the preconditioner.
  896.      * </p>
  897.      * <p>
  898.      * {@code shift} should be zero if the system A &middot; x = b is to be
  899.      * solved. Otherwise, it could be an approximation to an eigenvalue of A,
  900.      * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
  901.      * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
  902.      * sufficiently like an eigenvector corresponding to an eigenvalue near
  903.      * shift, then the computed x may have very large components. When
  904.      * normalized, x may be closer to an eigenvector than b.
  905.      * </p>
  906.      *
  907.      * @param a the linear operator A of the system
  908.      * @param m the preconditioner, M (can be {@code null})
  909.      * @param b the right-hand side vector
  910.      * @param goodb usually {@code false}, except if {@code x} is expected to
  911.      * contain a large multiple of {@code b}
  912.      * @param shift the amount to be subtracted to all diagonal elements of A
  913.      * @return a reference to {@code x} (shallow copy)
  914.      * @throws NullArgumentException if one of the parameters is {@code null}
  915.      * @throws MathIllegalArgumentException if {@code a} or {@code m} is not square
  916.      * @throws MathIllegalArgumentException if {@code m} or {@code b} have dimensions
  917.      * inconsistent with {@code a}
  918.      * @throws MathIllegalStateException at exhaustion of the iteration count,
  919.      * unless a custom
  920.      * {@link org.hipparchus.util.Incrementor.MaxCountExceededCallback callback}
  921.      * has been set at construction of the {@link IterationManager}
  922.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  923.      * {@code true}, and {@code a} or {@code m} is not self-adjoint
  924.      * @throws MathIllegalArgumentException if {@code m} is not
  925.      * positive definite
  926.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  927.      */
  928.     public RealVector solve(final RealLinearOperator a, final RealLinearOperator m,
  929.                             final RealVector b, final boolean goodb, final double shift)
  930.         throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException {
  931.         MathUtils.checkNotNull(a);
  932.         final RealVector x = new ArrayRealVector(a.getColumnDimension());
  933.         return solveInPlace(a, m, b, x, goodb, shift);
  934.     }

  936.     /**
  937.      * {@inheritDoc}
  938.      *
  939.      * @param x not meaningful in this implementation; should not be considered
  940.      * as an initial guess (<a href="#initguess">more</a>)
  941.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  942.      * {@code true}, and {@code a} or {@code m} is not self-adjoint
  943.      * @throws MathIllegalArgumentException if {@code m} is not positive
  944.      * definite
  945.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  946.      */
  947.     @Override
  948.     public RealVector solve(final RealLinearOperator a,
  949.         final RealLinearOperator m, final RealVector b, final RealVector x)
  950.         throws NullArgumentException,
  951.         MathIllegalArgumentException,
  952.         MathIllegalStateException {
  953.         MathUtils.checkNotNull(x);
  954.         return solveInPlace(a, m, b, x.copy(), false, 0.);
  955.     }

  957.     /**
  958.      * {@inheritDoc}
  959.      *
  960.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  961.      * {@code true}, and {@code a} is not self-adjoint
  962.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  963.      */
  964.     @Override
  965.     public RealVector solve(final RealLinearOperator a, final RealVector b)
  966.         throws NullArgumentException,
  967.         MathIllegalArgumentException, MathIllegalStateException {
  968.         MathUtils.checkNotNull(a);
  969.         final RealVector x = new ArrayRealVector(a.getColumnDimension());
  970.         x.set(0.);
  971.         return solveInPlace(a, null, b, x, false, 0.);
  972.     }

  974.     /**
  975.      * Returns the solution to the system (A - shift &middot; I) &middot; x = b.
  976.      * <p>
  977.      * If the solution x is expected to contain a large multiple of {@code b}
  978.      * (as in Rayleigh-quotient iteration), then better precision may be
  979.      * achieved with {@code goodb} set to {@code true}.
  980.      * </p>
  981.      * <p>
  982.      * {@code shift} should be zero if the system A &middot; x = b is to be
  983.      * solved. Otherwise, it could be an approximation to an eigenvalue of A,
  984.      * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
  985.      * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
  986.      * sufficiently like an eigenvector corresponding to an eigenvalue near
  987.      * shift, then the computed x may have very large components. When
  988.      * normalized, x may be closer to an eigenvector than b.
  989.      * </p>
  990.      *
  991.      * @param a the linear operator A of the system
  992.      * @param b the right-hand side vector
  993.      * @param goodb usually {@code false}, except if {@code x} is expected to
  994.      * contain a large multiple of {@code b}
  995.      * @param shift the amount to be subtracted to all diagonal elements of A
  996.      * @return a reference to {@code x}
  997.      * @throws NullArgumentException if one of the parameters is {@code null}
  998.      * @throws MathIllegalArgumentException if {@code a} is not square
  999.      * @throws MathIllegalArgumentException if {@code b} has dimensions
  1000.      * inconsistent with {@code a}
  1001.      * @throws MathIllegalStateException at exhaustion of the iteration count,
  1002.      * unless a custom
  1003.      * {@link org.hipparchus.util.Incrementor.MaxCountExceededCallback callback}
  1004.      * has been set at construction of the {@link IterationManager}
  1005.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  1006.      * {@code true}, and {@code a} is not self-adjoint
  1007.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  1008.      */
  1009.     public RealVector solve(final RealLinearOperator a, final RealVector b,
  1010.                             final boolean goodb, final double shift)
  1011.         throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException {
  1012.         MathUtils.checkNotNull(a);
  1013.         final RealVector x = new ArrayRealVector(a.getColumnDimension());
  1014.         return solveInPlace(a, null, b, x, goodb, shift);
  1015.     }

  1017.     /**
  1018.      * {@inheritDoc}
  1019.      *
  1020.      * @param x not meaningful in this implementation; should not be considered
  1021.      * as an initial guess (<a href="#initguess">more</a>)
  1022.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  1023.      * {@code true}, and {@code a} is not self-adjoint
  1024.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  1025.      */
  1026.     @Override
  1027.     public RealVector solve(final RealLinearOperator a, final RealVector b,
  1028.         final RealVector x) throws NullArgumentException, MathIllegalArgumentException,
  1029.         MathIllegalStateException {
  1030.         MathUtils.checkNotNull(x);
  1031.         return solveInPlace(a, null, b, x.copy(), false, 0.);
  1032.     }

  1034.     /**
  1035.      * {@inheritDoc}
  1036.      *
  1037.      * @param x the vector to be updated with the solution; {@code x} should
  1038.      * not be considered as an initial guess (<a href="#initguess">more</a>)
  1039.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  1040.      * {@code true}, and {@code a} or {@code m} is not self-adjoint
  1041.      * @throws MathIllegalArgumentException if {@code m} is not
  1042.      * positive definite
  1043.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  1044.      */
  1045.     @Override
  1046.     public RealVector solveInPlace(final RealLinearOperator a,
  1047.         final RealLinearOperator m, final RealVector b, final RealVector x)
  1048.         throws NullArgumentException,
  1049.         MathIllegalArgumentException,
  1050.         MathIllegalStateException {
  1051.         return solveInPlace(a, m, b, x, false, 0.);
  1052.     }

  1054.     /**
  1055.      * Returns an estimate of the solution to the linear system (A - shift
  1056.      * &middot; I) &middot; x = b. The solution is computed in-place.
  1057.      * <p>
  1058.      * If the solution x is expected to contain a large multiple of {@code b}
  1059.      * (as in Rayleigh-quotient iteration), then better precision may be
  1060.      * achieved with {@code goodb} set to {@code true}; this however requires an
  1061.      * extra call to the preconditioner.
  1062.      * </p>
  1063.      * <p>
  1064.      * {@code shift} should be zero if the system A &middot; x = b is to be
  1065.      * solved. Otherwise, it could be an approximation to an eigenvalue of A,
  1066.      * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
  1067.      * (b<sup>T</sup> &middot; b) corresponding to the vector b. If b is
  1068.      * sufficiently like an eigenvector corresponding to an eigenvalue near
  1069.      * shift, then the computed x may have very large components. When
  1070.      * normalized, x may be closer to an eigenvector than b.
  1071.      * </p>
  1072.      *
  1073.      * @param a the linear operator A of the system
  1074.      * @param m the preconditioner, M (can be {@code null})
  1075.      * @param b the right-hand side vector
  1076.      * @param x the vector to be updated with the solution; {@code x} should
  1077.      * not be considered as an initial guess (<a href="#initguess">more</a>)
  1078.      * @param goodb usually {@code false}, except if {@code x} is expected to
  1079.      * contain a large multiple of {@code b}
  1080.      * @param shift the amount to be subtracted to all diagonal elements of A
  1081.      * @return a reference to {@code x} (shallow copy).
  1082.      * @throws NullArgumentException if one of the parameters is {@code null}
  1083.      * @throws MathIllegalArgumentException if {@code a} or {@code m} is not square
  1084.      * @throws MathIllegalArgumentException if {@code m}, {@code b} or {@code x}
  1085.      * have dimensions inconsistent with {@code a}.
  1086.      * @throws MathIllegalStateException at exhaustion of the iteration count,
  1087.      * unless a custom
  1088.      * {@link org.hipparchus.util.Incrementor.MaxCountExceededCallback callback}
  1089.      * has been set at construction of the {@link IterationManager}
  1090.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  1091.      * {@code true}, and {@code a} or {@code m} is not self-adjoint
  1092.      * @throws MathIllegalArgumentException if {@code m} is not positive definite
  1093.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  1094.      */
  1095.     public RealVector solveInPlace(final RealLinearOperator a,
  1096.                                    final RealLinearOperator m, final RealVector b,
  1097.                                    final RealVector x, final boolean goodb, final double shift)
  1098.         throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException {
  1099.         checkParameters(a, m, b, x);

  1101.         final IterationManager manager = getIterationManager();
  1102.         /* Initialization counts as an iteration. */
  1103.         manager.resetIterationCount();
  1104.         manager.incrementIterationCount();

  1106.         final State state;
  1107.         state = new State(a, m, b, goodb, shift, delta, check);
  1108.         state.init();
  1109.         state.refineSolution(x);
  1110.         IterativeLinearSolverEvent event;
  1111.         event = new DefaultIterativeLinearSolverEvent(this,
  1112.                                                       manager.getIterations(),
  1113.                                                       x,
  1114.                                                       b,
  1115.                                                       state.getNormOfResidual());
  1116.         if (state.bEqualsNullVector()) {
  1117.             /* If b = 0 exactly, stop with x = 0. */
  1118.             manager.fireTerminationEvent(event);
  1119.             return x;
  1120.         }
  1121.         /* Cause termination if beta is essentially zero. */
  1122.         final boolean earlyStop;
  1123.         earlyStop = state.betaEqualsZero() || state.hasConverged();
  1124.         manager.fireInitializationEvent(event);
  1125.         if (!earlyStop) {
  1126.             do {
  1127.                 manager.incrementIterationCount();
  1128.                 event = new DefaultIterativeLinearSolverEvent(this,
  1129.                                                               manager.getIterations(),
  1130.                                                               x,
  1131.                                                               b,
  1132.                                                               state.getNormOfResidual());
  1133.                 manager.fireIterationStartedEvent(event);
  1134.                 state.update();
  1135.                 state.refineSolution(x);
  1136.                 event = new DefaultIterativeLinearSolverEvent(this,
  1137.                                                               manager.getIterations(),
  1138.                                                               x,
  1139.                                                               b,
  1140.                                                               state.getNormOfResidual());
  1141.                 manager.fireIterationPerformedEvent(event);
  1142.             } while (!state.hasConverged());
  1143.         }
  1144.         event = new DefaultIterativeLinearSolverEvent(this,
  1145.                                                       manager.getIterations(),
  1146.                                                       x,
  1147.                                                       b,
  1148.                                                       state.getNormOfResidual());
  1149.         manager.fireTerminationEvent(event);
  1150.         return x;
  1151.     }

  1153.     /**
  1154.      * {@inheritDoc}
  1155.      *
  1156.      * @param x the vector to be updated with the solution; {@code x} should
  1157.      * not be considered as an initial guess (<a href="#initguess">more</a>)
  1158.      * @throws MathIllegalArgumentException if {@link #shouldCheck()} is
  1159.      * {@code true}, and {@code a} is not self-adjoint
  1160.      * @throws MathIllegalArgumentException if {@code a} is ill-conditioned
  1161.      */
  1162.     @Override
  1163.     public RealVector solveInPlace(final RealLinearOperator a,
  1164.         final RealVector b, final RealVector x) throws NullArgumentException, MathIllegalArgumentException,
  1165.         MathIllegalStateException {
  1166.         return solveInPlace(a, null, b, x, false, 0.);
  1167.     }
  1168. }
Created with JaCoCo 0.8.12.202403310830