View Javadoc
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   */
17  
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;
23  
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;
31  
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 {
153 
154     /** {@code true} if symmetry of matrix and conditioner must be checked. */
155     private final boolean check;
156 
157     /**
158      * The value of the custom tolerance &delta; for the default stopping
159      * criterion.
160      */
161     private final double delta;
162 
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      */
231 
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;
251 
252         /** The machine precision. */
253         static final double MACH_PREC;
254 
255         /** Reference to the linear operator. */
256         private final RealLinearOperator a;
257 
258         /** Reference to the right-hand side vector. */
259         private final RealVector b;
260 
261         /** {@code true} if symmetry of matrix and conditioner must be checked. */
262         private final boolean check;
263 
264         /**
265          * The value of the custom tolerance &delta; for the default stopping
266          * criterion.
267          */
268         private final double delta;
269 
270         /** The value of beta[k+1]. */
271         private double beta;
272 
273         /** The value of beta[1]. */
274         private double beta1;
275 
276         /** The value of bstep[k-1]. */
277         private double bstep;
278 
279         /** The estimate of the norm of P * rC[k]. */
280         private double cgnorm;
281 
282         /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */
283         private double dbar;
284 
285         /**
286          * The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the
287          * initial code.
288          */
289         private double gammaZeta;
290 
291         /** The value of gbar[k]. */
292         private double gbar;
293 
294         /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
295         private double gmax;
296 
297         /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */
298         private double gmin;
299 
300         /** Copy of the {@code goodb} parameter. */
301         private final boolean goodb;
302 
303         /** {@code true} if the default convergence criterion is verified. */
304         private boolean hasConverged;
305 
306         /** The estimate of the norm of P * rL[k-1]. */
307         private double lqnorm;
308 
309         /** Reference to the preconditioner, M. */
310         private final RealLinearOperator m;
311 
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;
317 
318         /** The value of M * b. */
319         private final RealVector mb;
320 
321         /** The value of beta[k]. */
322         private double oldb;
323 
324         /** The value of beta[k] * M^(-1) * P' * v[k]. */
325         private RealVector r1;
326 
327         /** The value of beta[k+1] * M^(-1) * P' * v[k+1]. */
328         private RealVector r2;
329 
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;
335 
336         /** Copy of the {@code shift} parameter. */
337         private final double shift;
338 
339         /** The value of s[1] * ... * s[k-1]. */
340         private double snprod;
341 
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;
347 
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;
354 
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;
361 
362         /** The value of beta[k+1] * P' * v[k+1]. */
363         private RealVector y;
364 
365         /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */
366         private double ynorm2;
367 
368         /** The value of {@code b == 0} (exact floating-point equality). */
369         private boolean bIsNull;
370 
371         static {
372             MACH_PREC = FastMath.ulp(1.);
373             CBRT_MACH_PREC = FastMath.cbrt(MACH_PREC);
374         }
375 
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         }
408 
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         }
430 
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         }
439 
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         }
456 
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         }
477 
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         }
526 
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             }
543 
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;
609 
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         }
618 
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         }
743 
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);
759 
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         }
785 
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         }
794 
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         }
803 
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         }
813 
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     }
823 
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     }
840 
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     }
858 
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     }
869 
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         MathIllegalArgumentException, 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     }
887 
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     }
935 
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 MathIllegalArgumentException, NullArgumentException,
951         MathIllegalArgumentException,
952         MathIllegalStateException {
953         MathUtils.checkNotNull(x);
954         return solveInPlace(a, m, b, x.copy(), false, 0.);
955     }
956 
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 MathIllegalArgumentException, 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     }
973 
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     }
1016 
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 MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException,
1029         MathIllegalStateException {
1030         MathUtils.checkNotNull(x);
1031         return solveInPlace(a, null, b, x.copy(), false, 0.);
1032     }
1033 
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 MathIllegalArgumentException, NullArgumentException,
1049         MathIllegalArgumentException,
1050         MathIllegalStateException {
1051         return solveInPlace(a, m, b, x, false, 0.);
1052     }
1053 
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);
1100 
1101         final IterationManager manager = getIterationManager();
1102         /* Initialization counts as an iteration. */
1103         manager.resetIterationCount();
1104         manager.incrementIterationCount();
1105 
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     }
1152 
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 MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException,
1165         MathIllegalStateException {
1166         return solveInPlace(a, null, b, x, false, 0.);
1167     }
1168 }