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