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