Return-Path: X-Original-To: apmail-commons-commits-archive@minotaur.apache.org Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 618457F13 for ; Sat, 22 Oct 2011 13:01:41 +0000 (UTC) Received: (qmail 11607 invoked by uid 500); 22 Oct 2011 13:01:41 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 11372 invoked by uid 500); 22 Oct 2011 13:01:40 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 11365 invoked by uid 99); 22 Oct 2011 13:01:40 -0000 Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230) by apache.org (qpsmtpd/0.29) with ESMTP; Sat, 22 Oct 2011 13:01:40 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=5.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Sat, 22 Oct 2011 13:01:30 +0000 Received: from eris.apache.org (localhost [127.0.0.1]) by eris.apache.org (Postfix) with ESMTP id 049412388AA9 for ; Sat, 22 Oct 2011 13:01:08 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r1187709 [2/3] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/SymmLQ.java test/java/org/apache/commons/math/linear/SymmLQTest.java Date: Sat, 22 Oct 2011 13:01:07 -0000 To: commits@commons.apache.org From: celestin@apache.org X-Mailer: svnmailer-1.0.8-patched Message-Id: <20111022130108.049412388AA9@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java?rev=1187709&r1=1187708&r2=1187709&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java Sat Oct 22 13:01:07 2011 @@ -1,1223 +1,1223 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math.linear; - -import org.apache.commons.math.exception.DimensionMismatchException; -import org.apache.commons.math.exception.MaxCountExceededException; -import org.apache.commons.math.exception.NullArgumentException; -import org.apache.commons.math.exception.util.ExceptionContext; -import org.apache.commons.math.util.FastMath; -import org.apache.commons.math.util.IterationManager; -import org.apache.commons.math.util.MathUtils; - -/** - *

- * Implementation of the SYMMLQ iterative linear solver proposed by Paige and Saunders (1975). This implementation is - * largely based on the FORTRAN code by Pr. Michael A. Saunders, available here. - *

- *

- * SYMMLQ is designed to solve the system of linear equations A · x = b - * where A is an n × n self-adjoint linear operator (defined as a - * {@link RealLinearOperator}), and b is a given vector. The operator A is not - * required to be positive definite. If A is known to be definite, the method of - * conjugate gradients might be preferred, since it will require about the same - * number of iterations as SYMMLQ but slightly less work per iteration. - *

- *

- * SYMMLQ is designed to solve the system (A - shift · I) · x = b, - * where shift is a specified scalar value. If shift and b are suitably chosen, - * the computed vector x may approximate an (unnormalized) eigenvector of A, as - * in the methods of inverse iteration and/or Rayleigh-quotient iteration. - * Again, the linear operator (A - shift · I) need not be positive - * definite (but must be self-adjoint). The work per iteration is very - * slightly less if shift = 0. - *

- *

Peconditioning

- *

- * Preconditioning may reduce the number of iterations required. The solver is - * provided with a positive definite preconditioner M = C · C^T - * that is known to approximate (A - shift · I) in some sense, while - * systems of the form M · y = x can be solved efficiently. Then SYMMLQ - * will implicitly solve the system of equations P · (A - shift · - * I) · P^T · xhat = P · b, i.e. Ahat · - * xhat = bhat, where P = C^-1, Ahat = P · (A - shift · - * I) · P^T, bhat = P · b, and return the solution x = - * P^T · xhat. The associated residual is rhat = bhat - Ahat - * · xhat = P · [b - (A - shift · I) · x] = P - * · r. - *

- *

Default stopping criterion

- *

- * A default stopping criterion is implemented. The iterations stop when || rhat - * || ≤ δ || Ahat || || xhat ||, where xhat is the current estimate of - * the solution of the transformed system, rhat the current estimate of the - * corresponding residual, and δ a user-specified tolerance. - *

- *

Iteration count

- *

- * In the present context, an iteration should be understood as one evaluation - * of the matrix-vector product A · x. The initialization phase therefore - * counts as one iteration. If the user requires checks on the symmetry of A, - * this entails one further matrix-vector product by iteration. This further - * product is not accounted for in the iteration count. In other words, - * the number of iterations required to reach convergence will be identical, - * whether checks have been required or not. - *

- *

- * The present definition of the iteration count differs from that adopted in - * the original FOTRAN code, where the initialization phase was not - * taken into account. - *

- *

Exception context

- *

- * Besides standard {@link DimensionMismatchException}, this class might throw - * {@link NonSelfAdjointOperatorException} if the linear operator or the - * preconditioner are not symmetric. In this case, the {@link ExceptionContext} - * provides more information - *

key {@code "operator"} points to the offending linear operator, say L,
key {@code "vector1"} points to the first offending vector, say x, - *
key {@code "vector2"} points to the second offending vector, say y, such - * that x^T · L · y ≠ y^T · L - * · x (within a certain accuracy).

- *

- * {@link NonPositiveDefiniteOperatorException} might also be thrown in case the - * preconditioner is not positive definite. The relevant keys to the - * {@link ExceptionContext} are - *

key {@code "operator"}, which points to the offending linear operator, - * say L,
key {@code "vector"}, which points to the offending vector, say x, such - * that x^T · L · x < 0.

- *

References

- *

Paige and Saunders (1975): C. C. Paige and M. A. Saunders, - * Solution of Sparse Indefinite Systems of Linear Equations, SIAM - * Journal on Numerical Analysis 12(4): 617-629, 1975

- * - * @version $Id$ - * @since 3.0 - */ -public class SymmLQ - extends PreconditionedIterativeLinearSolver { - - /* - * IMPLEMENTATION NOTES - * -------------------- - * The implementation follows as closely as possible the notations of Paige - * and Saunders (1975). Attention must be paid to the fact that some - * quantities which are relevant to iteration k can only be computed in - * iteration (k+1). Therefore, minute attention must be paid to the index of - * each state variable of this algorithm. - * - * 1. Preconditioning - * --------------- - * The Lanczos iterations associated with Ahat and bhat read - * beta[1] = |P . b| - * v[1] = P.b / beta[1] - * beta[k+1] * v[k+1] = Ahat * v[k] - alpha[k] * v[k] - beta[k] * v[k-1] - * = P * (A - shift * I) * P' * v[k] - alpha[k] * v[k] - * - beta[k] * v[k-1] - * Multiplying both sides by P', we get - * beta[k+1] * (P' * v)[k+1] = M^(-1) * (A - shift * I) * (P' * v)[k] - * - alpha[k] * (P' * v)[k] - * - beta[k] * (P' * v[k-1]), - * and - * alpha[k+1] = v[k+1]' * Ahat * v[k+1] - * = v[k+1]' * P * (A - shift * I) * P' * v[k+1] - * = (P' * v)[k+1]' * (A - shift * I) * (P' * v)[k+1]. - * - * In other words, the Lanczos iterations are unchanged, except for the fact - * that we really compute (P' * v) instead of v. It can easily be checked - * that all other formulas are unchanged. It must be noted that P is never - * explicitly used, only matrix-vector products involving M^(-1) are - * invoked. - * - * 2. Accounting for the shift parameter - * ---------------------------------- - * Is trivial: each time A.operate(x) is invoked, one must subtract shift * x - * to the result. - * - * 3. Accounting for the goodb flag - * ----------------------------- - * When goodb is set to true, the component of xL along b is computed - * separately. From Page and Saunders (1975), equation (5.9), we have - * wbar[k+1] = s[k] * wbar[k] - c[k] * v[k+1], - * wbar[1] = v[1]. - * Introducing wbar2[k] = wbar[k] - s[1] * ... * s[k-1] * v[1], it can - * easily be verified by induction that what follows the same recursive - * relation - * wbar2[k+1] = s[k] * wbar2[k] - c[k] * v[k+1], - * wbar2[1] = 0, - * and we then have - * w[k] = c[k] * wbar2[k] + s[k] * v[k+1] - * + s[1] * ... * s[k-1] * c[k] * v[1]. - * Introducing w2[k] = w[k] - s[1] * ... * s[k-1] * c[k] * v[1], we find, - * from (5.10) - * xL[k] = zeta[1] * w[1] + ... + zeta[k] * w[k] - * = zeta[1] * w2[1] + ... + zeta[k] * w2[k] - * + (s[1] * c[2] * zeta[2] + ... - * + s[1] * ... * s[k-1] * c[k] * zeta[k]) * v[1] - * = xL2[k] + bstep[k] * v[1], - * where xL2[k] is defined by - * xL2[0] = 0, - * xL2[k+1] = xL2[k] + zeta[k+1] * w2[k+1], - * and bstep is defined by - * bstep[1] = 0, - * bstep[k] = bstep[k-1] + s[1] * ... * s[k-1] * c[k] * zeta[k]. - * We also have, from (5.11) - * xC[k] = xL[k-1] + zbar[k] * wbar[k] - * = xL2[k-1] + zbar[k] * wbar2[k] - * + (bstep[k-1] + s[1] * ... * s[k-1] * zbar[k]) * v[1]. - */ - - /** - *

- * A simple container holding the non-final variables used in the - * iterations. Making the current state of the solver visible from the - * outside is necessary, because during the iterations, {@code x} does not - * exactly hold the current estimate of the solution. Indeed, - * {@code x} needs in general to be moved from the LQ point to the CG point. - * Besides, additional upudates must be carried out in case {@code goodb} is - * set to {@code true}. - *

- *

- * In all subsequent comments, the description of the state variables refer - * to their value after a call to {@link #update()}. In these comments, k is - * the current number of evaluations of matrix-vector products. - *

- */ - private class State { - - /** Reference to the linear operator. */ - private final RealLinearOperator a; - - /** Reference to the right-hand side vector. */ - private final RealVector b; - - /** The value of beta[k+1]. */ - private double beta; - - /** The value of beta[1]. */ - private double beta1; - - /** The value of bstep[k-1]. */ - private double bstep; - - /** The estimate of the norm of P * rC[k]. */ - private double cgnorm; - - /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */ - private double dbar; - - /** - * The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the - * initial code. - */ - private double gammaZeta; - - /** The value of gbar[k]. */ - private double gbar; - - /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */ - private double gmax; - - /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */ - private double gmin; - - /** Copy of the {@code goodb} parameter. */ - private final boolean goodb; - - /** {@code true} if the default convergence criterion is verified. */ - private boolean hasConverged; - - /** The estimate of the norm of P * rL[k-1]. */ - private double lqnorm; - - /** Reference to the preconditioner. */ - private final InvertibleRealLinearOperator m; - - /** - * The value of (-eps[k+1] * zeta[k-1]). Was called {@code rhs2} in the - * initial code. - */ - private double minusEpsZeta; - - /** The value of M^(-1) * b. */ - private final RealVector mSolveB; - - /** The value of beta[k]. */ - private double oldb; - - /** The value of beta[k] * M * P' * v[k]. */ - private RealVector r1; - - /** The value of beta[k+1] * M * P' * v[k+1]. */ - private RealVector r2; - - /** Copy of the {@code shift} parameter. */ - private final double shift; - - /** The value of s[1] * ... * s[k-1]. */ - private double snprod; - - /** - * An estimate of the square of the norm of A * V[k], based on Paige and - * Saunders (1975), equation (3.3). - */ - private double tnorm; - - /** - * The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] * - * v[1]) if {@code goodb} is {@code true}. Was called {@code w} in the - * initial code. - */ - private RealVector wbar; - - /** - * A reference to the vector to be updated with the solution. Contains - * the value of xL[k-1] if {@code goodb} is {@code false}, (xL[k-1] - - * bstep[k-1] * v[1]) otherwise. - */ - private final RealVector x; - - /** The value of beta[k+1] * P' * v[k+1]. */ - private RealVector y; - - /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */ - private double ynorm2; - - /** - * Creates and inits to k = 1 a new instance of this class. - * - * @param a Linear operator A of the system. - * @param m Preconditioner (can be {@code null}). - * @param b Right-hand side vector. - * @param x Vector to be updated with the solution. {@code x} should not - * be considered as an initial guess, as it is set to 0 in the - * initialization phase. - * @param goodb Usually {@code false}, except if {@code x} is expected - * to contain a large multiple of {@code b}. - * @param shift The amount to be subtracted to all diagonal elements of - * A. - */ - public State(final RealLinearOperator a, - final InvertibleRealLinearOperator m, final RealVector b, - final RealVector x, final boolean goodb, final double shift) { - this.a = a; - this.m = m; - this.b = b; - this.x = x; - this.goodb = goodb; - this.shift = shift; - this.mSolveB = m == null ? b : m.solve(b); - this.hasConverged = false; - init(); - } - - /** - * Move to the CG point if it seems better. In this version of SYMMLQ, - * the convergence tests involve only cgnorm, so we're unlikely to stop - * at an LQ point, except if the iteration limit interferes. - * - * @param xRefined Vector to be updated with the refined value of x. - */ - public void refine(final RealVector xRefined) { - final int n = this.x.getDimension(); - if (lqnorm < cgnorm) { - if (!goodb) { - xRefined.setSubVector(0, this.x); - } else { - final double step = bstep / beta1; - for (int i = 0; i < n; i++) { - final double bi = mSolveB.getEntry(i); - final double xi = this.x.getEntry(i); - xRefined.setEntry(i, xi + step * bi); - } - } - } else { - final double anorm = FastMath.sqrt(tnorm); - final double diag = gbar == 0. ? anorm * MACH_PREC : gbar; - final double zbar = gammaZeta / diag; - final double step = (bstep + snprod * zbar) / beta1; - // ynorm = FastMath.sqrt(ynorm2 + zbar * zbar); - if (!goodb) { - for (int i = 0; i < n; i++) { - final double xi = this.x.getEntry(i); - final double wi = wbar.getEntry(i); - xRefined.setEntry(i, xi + zbar * wi); - } - } else { - for (int i = 0; i < n; i++) { - final double xi = this.x.getEntry(i); - final double wi = wbar.getEntry(i); - final double bi = mSolveB.getEntry(i); - xRefined.setEntry(i, xi + zbar * wi + step * bi); - } - } - } - } - - /** - * Performs the initial phase of the SYMMLQ algorithm. On return, the - * value of the state variables of {@code this} object correspond to k = - * 1. - */ - private void init() { - this.x.set(0.); - /* - * Set up y for the first Lanczos vector. y and beta1 will be zero - * if b = 0. - */ - this.r1 = this.b.copy(); - this.y = this.m == null ? this.b.copy() : this.m.solve(this.r1); - if ((this.m != null) && check) { - checkSymmetry(this.m, this.r1, this.y, this.m.solve(this.y)); - } - - this.beta1 = this.r1.dotProduct(this.y); - if (this.beta1 < 0.) { - throwNPDLOException(this.m, this.y); - } - if (this.beta1 == 0.) { - /* If b = 0 exactly, stop with x = 0. */ - return; - } - this.beta1 = FastMath.sqrt(this.beta1); - /* At this point - * r1 = b, - * y = M^(-1) * b, - * beta1 = beta[1]. - */ - final RealVector v = this.y.mapMultiply(1. / this.beta1); - this.y = this.a.operate(v); - if (check) { - checkSymmetry(this.a, v, this.y, this.a.operate(this.y)); - } - /* - * Set up y for the second Lanczos vector. y and beta will be zero - * or very small if b is an eigenvector. - */ - daxpy(-this.shift, v, this.y); - final double alpha = v.dotProduct(this.y); - daxpy(-alpha / this.beta1, this.r1, this.y); - /* - * At this point - * alpha = alpha[1] - * y = beta[2] * M * P' * v[2] - */ - /* Make sure r2 will be orthogonal to the first v. */ - final double vty = v.dotProduct(this.y); - final double vtv = v.dotProduct(v); - daxpy(-vty / vtv, v, this.y); - this.r2 = this.y.copy(); - if (this.m != null) { - this.y = this.m.solve(this.r2); - } - this.oldb = this.beta1; - this.beta = this.r2.dotProduct(this.y); - if (this.beta < 0.) { - throwNPDLOException(this.m, this.y); - } - this.beta = FastMath.sqrt(this.beta); - /* - * At this point - * oldb = beta[1] - * beta = beta[2] - * y = beta[2] * P' * v[2] - * r2 = beta[2] * M * P' * v[2] - */ - this.cgnorm = this.beta1; - this.gbar = alpha; - this.dbar = this.beta; - this.gammaZeta = this.beta1; - this.minusEpsZeta = 0.; - this.bstep = 0.; - this.snprod = 1.; - this.tnorm = alpha * alpha + this.beta * this.beta; - this.ynorm2 = 0.; - this.gmax = FastMath.abs(alpha) + MACH_PREC; - this.gmin = this.gmax; - - if (this.goodb) { - this.wbar = new ArrayRealVector(this.a.getRowDimension()); - this.wbar.set(0.); - } else { - this.wbar = v; - } - updateNorms(); - } - - /** - * Performs the next iteration of the algorithm. The iteration count - * should be incremented prior to calling this method. On return, the - * value of the state variables of {@code this} object correspond to the - * current iteration count {@code k}. - */ - private void update() { - final RealVector v = y.mapMultiply(1. / beta); - y = a.operate(v); - daxpbypz(-shift, v, -beta / oldb, r1, y); - final double alpha = v.dotProduct(y); - /* - * At this point - * v = P' * v[k], - * y = (A - shift * I) * P' * v[k] - beta[k] * M * P' * v[k-1], - * alpha = v'[k] * P * (A - shift * I) * P' * v[k] - * - beta[k] * v[k]' * P * M * P' * v[k-1] - * = v'[k] * P * (A - shift * I) * P' * v[k] - * - beta[k] * v[k]' * v[k-1] - * = alpha[k]. - */ - daxpy(-alpha / beta, r2, y); - /* - * At this point - * y = (A - shift * I) * P' * v[k] - alpha[k] * M * P' * v[k] - * - beta[k] * M * P' * v[k-1] - * = M * P' * (P * (A - shift * I) * P' * v[k] -alpha[k] * v[k] - * - beta[k] * v[k-1]) - * = beta[k+1] * M * P' * v[k+1], - * from Paige and Saunders (1975), equation (3.2). - * - * WATCH-IT: the two following line work only because y is no longer - * updated up to the end of the present iteration, and is - * reinitialized at the beginning of the next iteration. - */ - r1 = r2; - r2 = y; - if (m != null) { - y = m.solve(r2); - } - oldb = beta; - beta = r2.dotProduct(y); - if (beta < 0.) { - throwNPDLOException(m, y); - } - beta = FastMath.sqrt(beta); - /* - * At this point - * r1 = beta[k] * M * P' * v[k], - * r2 = beta[k+1] * M * P' * v[k+1], - * y = beta[k+1] * P' * v[k+1], - * oldb = beta[k], - * beta = beta[k+1]. - */ - tnorm += alpha * alpha + oldb * oldb + beta * beta; - /* - * Compute the next plane rotation for Q. See Paige and Saunders - * (1975), equation (5.6), with - * gamma = gamma[k-1], - * c = c[k-1], - * s = s[k-1]. - */ - final double gamma = FastMath.sqrt(gbar * gbar + oldb * oldb); - final double c = gbar / gamma; - final double s = oldb / gamma; - /* - * The relations - * gbar[k] = s[k-1] * (-c[k-2] * beta[k]) - c[k-1] * alpha[k] - * = s[k-1] * dbar[k] - c[k-1] * alpha[k], - * delta[k] = c[k-1] * dbar[k] + s[k-1] * alpha[k], - * are not stated in Paige and Saunders (1975), but can be retrieved - * by expanding the (k, k-1) and (k, k) coefficients of the matrix in - * equation (5.5). - */ - final double deltak = c * dbar + s * alpha; - gbar = s * dbar - c * alpha; - final double eps = s * beta; - dbar = -c * beta; - final double zeta = gammaZeta / gamma; - /* - * At this point - * gbar = gbar[k] - * deltak = delta[k] - * eps = eps[k+1] - * dbar = dbar[k+1] - * zeta = zeta[k-1] - */ - final double zetaC = zeta * c; - final double zetaS = zeta * s; - final int n = x.getDimension(); - for (int i = 0; i < n; i++) { - final double xi = x.getEntry(i); - final double vi = v.getEntry(i); - final double wi = wbar.getEntry(i); - x.setEntry(i, xi + wi * zetaC + vi * zetaS); - wbar.setEntry(i, wi * s - vi * c); - } - /* - * At this point - * x = xL[k-1], - * ptwbar = P' wbar[k], - * see Paige and Saunders (1975), equations (5.9) and (5.10). - */ - bstep += snprod * c * zeta; - snprod *= s; - gmax = FastMath.max(gmax, gamma); - gmin = FastMath.min(gmin, gamma); - ynorm2 += zeta * zeta; - gammaZeta = minusEpsZeta - deltak * zeta; - minusEpsZeta = -eps * zeta; - /* - * At this point - * snprod = s[1] * ... * s[k-1], - * gmax = max(|alpha[1]|, gamma[1], ..., gamma[k-1]), - * gmin = min(|alpha[1]|, gamma[1], ..., gamma[k-1]), - * ynorm2 = zeta[1]^2 + ... + zeta[k-1]^2, - * gammaZeta = gamma[k] * zeta[k], - * minusEpsZeta = -eps[k+1] * zeta[k-1]. - * The relation for gammaZeta can be retrieved from Paige and - * Saunders (1975), equation (5.4a), last line of the vector - * gbar[k] * zbar[k] = -eps[k] * zeta[k-2] - delta[k] * zeta[k-1]. - */ - updateNorms(); - } - - /** - * Computes the norms of the residuals, and checks for convergence. - * Updates {@link #lqnorm} and {@link #cgnorm}. - */ - private void updateNorms() { - final double anorm = FastMath.sqrt(tnorm); - final double ynorm = FastMath.sqrt(ynorm2); - final double epsa = anorm * MACH_PREC; - final double epsx = anorm * ynorm * MACH_PREC; - final double epsr = anorm * ynorm * delta; - final double diag = gbar == 0. ? epsa : gbar; - lqnorm = FastMath.sqrt(gammaZeta * gammaZeta + - minusEpsZeta * minusEpsZeta); - final double qrnorm = snprod * beta1; - cgnorm = qrnorm * beta / FastMath.abs(diag); - - /* - * Estimate cond(A). In this version we look at the diagonals of L - * in the factorization of the tridiagonal matrix, T = L * Q. - * Sometimes, T[k] can be misleadingly ill-conditioned when T[k+1] - * is not, so we must be careful not to overestimate acond. - */ - final double acond; - if (lqnorm <= cgnorm) { - acond = gmax / gmin; - } else { - acond = gmax / FastMath.min(gmin, FastMath.abs(diag)); - } - if (acond * MACH_PREC >= 0.1) { - throw new IllConditionedOperatorException(acond); - } - if (beta1 <= epsx) { - /* - * x has converged to an eigenvector of A corresponding to the - * eigenvalue shift. - */ - throw new SingularOperatorException(); - } - hasConverged = (cgnorm <= epsx) || (cgnorm <= epsr); - } - } - - /** The cubic root of {@link #MACH_PREC}. */ - private static final double CBRT_MACH_PREC; - - /** The machine precision. */ - private static final double MACH_PREC; - - /** Key for the exception context. */ - private static final String OPERATOR = "operator"; - - /** Key for the exception context. */ - private static final String THRESHOLD = "threshold"; - - /** Key for the exception context. */ - private static final String VECTOR = "vector"; - - /** Key for the exception context. */ - private static final String VECTOR1 = "vector1"; - - /** Key for the exception context. */ - private static final String VECTOR2 = "vector2"; - - /** {@code true} if symmetry of matrix and conditioner must be checked. */ - private final boolean check; - - /** - * The value of the custom tolerance δ for the default stopping - * criterion. - */ - private final double delta; - - /** - * Creates a new instance of this class, with default - * stopping criterion. - * - * @param maxIterations Maximum number of iterations. - * @param delta δ parameter for the default stopping criterion. - * @param check {@code true} if self-adjointedness of both matrix and - * preconditioner should be checked. This entails an extra matrix-vector - * product at each iteration. - */ - public SymmLQ(final int maxIterations, final double delta, - final boolean check) { - super(maxIterations); - this.delta = delta; - this.check = check; - } - - /** - * Creates a new instance of this class, with default - * stopping criterion and custom iteration manager. - * - * @param manager Custom iteration manager. - * @param delta δ parameter for the default stopping criterion. - * @param check {@code true} if self-adjointedness of both matrix and - * preconditioner should be checked. This entails an extra matrix-vector - * product at each iteration. - */ - public SymmLQ(final IterationManager manager, final double delta, - final boolean check) { - super(manager); - this.delta = delta; - this.check = check; - } - - static { - MACH_PREC = Math.ulp(1.); - CBRT_MACH_PREC = Math.cbrt(MACH_PREC); - } - - /** - * Performs a symmetry check on the specified linear operator, and throws an - * exception in case this check fails. Given a linear operator L, and a - * vector x, this method checks that x' L y = y' L x (within a given - * accuracy), where y = L x. - * - * @param l The linear operator L. - * @param x The candidate vector x. - * @param y The candidate vector y = L x. - * @param z The vector z = L y. - * @throws NonSelfAdjointOperatorException when the test fails. - */ - private static void checkSymmetry(final RealLinearOperator l, - final RealVector x, final RealVector y, - final RealVector z) - throws NonSelfAdjointOperatorException { - final double s = y.dotProduct(y); - final double t = x.dotProduct(z); - final double epsa = (s + MACH_PREC) * CBRT_MACH_PREC; - if (FastMath.abs(s - t) > epsa) { - final NonSelfAdjointOperatorException e; - e = new NonSelfAdjointOperatorException(); - final ExceptionContext context = e.getContext(); - context.setValue(OPERATOR, l); - context.setValue(VECTOR1, x); - context.setValue(VECTOR2, y); - context.setValue(THRESHOLD, Double.valueOf(epsa)); - throw e; - } - } - - /** - * A BLAS-like function, for the operation z ← a · x + b - * · y + z. This is for internal use only: no dimension checks are - * provided. - * - * @param a The scalar by which {@code x} is to be multiplied. - * @param x The first vector to be added to {@code z}. - * @param b The scalar by which {@code y} is to be multiplied. - * @param y The second vector to be added to {@code z}. - * @param z The vector to be incremented. - */ - private static void daxpbypz(final double a, final RealVector x, - final double b, final RealVector y, - final RealVector z) { - final int n = z.getDimension(); - for (int i = 0; i < n; i++) { - final double zi; - zi = a * x.getEntry(i) + b * y.getEntry(i) + z.getEntry(i); - z.setEntry(i, zi); - } - } - - /** - * A clone of the BLAS {@code DAXPY} function, which carries out the - * operation y ← a · x + y. This is for internal use only: no - * dimension checks are provided. - * - * @param a The scalar by which {@code x} is to be multiplied. - * @param x The vector to be added to {@code y}. - * @param y The vector to be incremented. - */ - private static void daxpy(final double a, final RealVector x, - final RealVector y) { - final int n = x.getDimension(); - for (int i = 0; i < n; i++) { - y.setEntry(i, a * x.getEntry(i) + y.getEntry(i)); - } - } - - /** - * Throws a new {@link NonPositiveDefiniteOperatorException} with - * appropriate context. - * - * @param l The offending linear operator. - * @param v The offending vector. - * @throws NonPositiveDefiniteOperatorException in any circumstances. - */ - private static void throwNPDLOException(final RealLinearOperator l, - final RealVector v) - throws NonPositiveDefiniteOperatorException { - final NonPositiveDefiniteOperatorException e; - e = new NonPositiveDefiniteOperatorException(); - final ExceptionContext context = e.getContext(); - context.setValue(OPERATOR, l); - context.setValue(VECTOR, v); - throw e; - } - - /** - * Returns {@code true} if symmetry of the matrix, and symmetry as well as - * positive definiteness of the preconditioner should be checked. - * - * @return {@code true} if the tests are to be performed. - */ - public final boolean getCheck() { - return check; - } - - /** - * Returns an estimate of the solution to the linear system A · x = - * b. - * - * @param a Linear operator A of the system. - * @param m Preconditioner (can be {@code null}). - * @param b Right-hand side vector. - * @return A new vector containing the solution. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} or {@code m} is not - * square. - * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x} - * have dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} or {@code m} is not self-adjoint. - * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive - * definite. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - @Override - public RealVector solve(final RealLinearOperator a, - final InvertibleRealLinearOperator m, - final RealVector b) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - NonPositiveDefiniteOperatorException, IllConditionedOperatorException, - MaxCountExceededException { - MathUtils.checkNotNull(a); - final RealVector x = new ArrayRealVector(a.getColumnDimension()); - return solveInPlace(a, m, b, x, false, 0.); - } - - /** - * Returns an estimate of the solution to the linear system (A - shift - * · I) · x = b. - *

- * If the solution x is expected to contain a large multiple of {@code b} - * (as in Rayleigh-quotient iteration), then better precision may be - * achieved with {@code goodb} set to {@code true}; this however requires an - * extra call to the preconditioner. - *

- *

- * {@code shift} should be zero if the system A · x = b is to be - * solved. Otherwise, it could be an approximation to an eigenvalue of A, - * such as the Rayleigh quotient b^T · A · b / - * (b^T · b) corresponding to the vector b. If b is - * sufficiently like an eigenvector corresponding to an eigenvalue near - * shift, then the computed x may have very large components. When - * normalized, x may be closer to an eigenvector than b. - *

- * - * @param a Linear operator A of the system. - * @param m Preconditioner (can be {@code null}). - * @param b Right-hand side vector. - * @param goodb Usually {@code false}, except if {@code x} is expected to - * contain a large multiple of {@code b}. - * @param shift The amount to be subtracted to all diagonal elements of A. - * @return A reference to {@code x} (shallow copy). - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} or {@code m} is not - * square. - * @throws DimensionMismatchException if {@code m} or {@code b} have - * dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} or {@code m} is not self-adjoint. - * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive - * definite. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - public RealVector solve(final RealLinearOperator a, - final InvertibleRealLinearOperator m, - final RealVector b, final boolean goodb, - final double shift) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - NonPositiveDefiniteOperatorException, IllConditionedOperatorException, - MaxCountExceededException { - MathUtils.checkNotNull(a); - final RealVector x = new ArrayRealVector(a.getColumnDimension()); - return solveInPlace(a, m, b, x, goodb, shift); - } - - /** - * Returns an estimate of the solution to the linear system A · x = - * b. - * - * @param a Linear operator A of the system. - * @param m Preconditioner (can be {@code null}). - * @param b Right-hand side vector. - * @param x Not meaningful in this implementation. Should not be considered - * as an initial guess. - * @return A new vector containing the solution. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} or {@code m} is not - * square. - * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x} - * have dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} or {@code m} is not self-adjoint. - * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive - * definite. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - @Override - public RealVector solve(final RealLinearOperator a, - final InvertibleRealLinearOperator m, - final RealVector b, final RealVector x) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - NonPositiveDefiniteOperatorException, IllConditionedOperatorException, - MaxCountExceededException { - MathUtils.checkNotNull(x); - return solveInPlace(a, m, b, x.copy(), false, 0.); - } - - /** - * Returns an estimate of the solution to the linear system A · x = - * b. - * - * @param a Linear operator A of the system. - * @param b Right-hand side vector. - * @return A new vector containing the solution. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} is not square. - * @throws DimensionMismatchException if {@code b} has dimensions - * inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} is not self-adjoint. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - @Override - public RealVector solve(final RealLinearOperator a, final RealVector b) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - IllConditionedOperatorException, MaxCountExceededException { - MathUtils.checkNotNull(a); - final RealVector x = new ArrayRealVector(a.getColumnDimension()); - x.set(0.); - return solveInPlace(a, null, b, x, false, 0.); - } - - /** - * Returns the solution to the system (A - shift · I) · x = b. - *

- *

- * - * @param a Linear operator A of the system. - * @param b Right-hand side vector. - * @param goodb Usually {@code false}, except if {@code x} is expected to - * contain a large multiple of {@code b}. - * @param shift The amount to be subtracted to all diagonal elements of A. - * @return a reference to {@code x}. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} is not square. - * @throws DimensionMismatchException if {@code b} has dimensions - * inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} is not self-adjoint. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - public RealVector solve(final RealLinearOperator a, final RealVector b, - final boolean goodb, final double shift) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - IllConditionedOperatorException, MaxCountExceededException { - MathUtils.checkNotNull(a); - final RealVector x = new ArrayRealVector(a.getColumnDimension()); - return solveInPlace(a, null, b, x, goodb, shift); - } - - /** - * Returns an estimate of the solution to the linear system A · x = - * b. - * - * @param a Linear operator A of the system. - * @param b Right-hand side vector. - * @param x Not meaningful in this implementation. Should not be considered - * as an initial guess. - * @return A new vector containing the solution. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} is not square. - * @throws DimensionMismatchException if {@code b} or {@code x} have - * dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} is not self-adjoint. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - @Override - public RealVector solve(final RealLinearOperator a, final RealVector b, - final RealVector x) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - IllConditionedOperatorException, MaxCountExceededException { - MathUtils.checkNotNull(x); - return solveInPlace(a, null, b, x.copy(), false, 0.); - } - - /** - * Returns an estimate of the solution to the linear system A · x = - * b. The solution is computed in-place. - * - * @param a Linear operator A of the system. - * @param m Preconditioner (can be {@code null}). - * @param b Right-hand side vector. - * @param x Vector to be updated with the solution. {@code x} should not be - * considered as an initial guess, as it is set to 0 in the initialization - * phase. - * @return A reference to {@code x} (shallow copy) updated with the - * solution. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} or {@code m} is not - * square. - * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x} - * have dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} or {@code m} is not self-adjoint. - * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive - * definite. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - @Override - public RealVector solveInPlace(final RealLinearOperator a, - final InvertibleRealLinearOperator m, - final RealVector b, final RealVector x) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - NonPositiveDefiniteOperatorException, IllConditionedOperatorException, - MaxCountExceededException { - return solveInPlace(a, m, b, x, false, 0.); - } - - /** - * Returns an estimate of the solution to the linear system (A - shift - * · I) · x = b. The solution is computed in-place. - *

- *

- * - * @param a Linear operator A of the system. - * @param m Preconditioner (can be {@code null}). - * @param b Right-hand side vector. - * @param x Vector to be updated with the solution. {@code x} should not be - * considered as an initial guess, as it is set to 0 in the initialization - * phase. - * @param goodb Usually {@code false}, except if {@code x} is expected to - * contain a large multiple of {@code b}. - * @param shift The amount to be subtracted to all diagonal elements of A. - * @return A reference to {@code x} (shallow copy). - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} or {@code m} is not - * square. - * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x} - * have dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} or {@code m} is not self-adjoint. - * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive - * definite. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - public RealVector solveInPlace(final RealLinearOperator a, - final InvertibleRealLinearOperator m, - final RealVector b, final RealVector x, - final boolean goodb, final double shift) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - NonPositiveDefiniteOperatorException, IllConditionedOperatorException, - MaxCountExceededException { - checkParameters(a, m, b, x); - - final IterationManager manager = getIterationManager(); - /* Initialization counts as an iteration. */ - manager.resetIterationCount(); - manager.incrementIterationCount(); - - final State state = new State(a, m, b, x, goodb, shift); - final IterativeLinearSolverEvent event = createEvent(state); - if (state.beta1 == 0.) { - /* If b = 0 exactly, stop with x = 0. */ - manager.fireTerminationEvent(event); - return x; - } - /* Cause termination if beta is essentially zero. */ - final boolean earlyStop; - earlyStop = (state.beta < MACH_PREC) || (state.hasConverged); - manager.fireInitializationEvent(event); - if (!earlyStop) { - do { - manager.incrementIterationCount(); - manager.fireIterationStartedEvent(event); - state.update(); - manager.fireIterationPerformedEvent(event); - } while (!state.hasConverged); - } - state.refine(x); - /* - * The following two lines are a hack because state.x is now refined, - * so further calls to state.refine() (via event.getSolution()) should - * *not* return an altered value of state.x. - */ - state.bstep = 0.; - state.gammaZeta = 0.; - manager.fireTerminationEvent(event); - return x; - } - - /** - * Returns an estimate of the solution to the linear system A · x = - * b. The solution is computed in-place. - * - * @param a Linear operator A of the system. - * @param b Right-hand side vector. - * @param x Vector to be updated with the solution. {@code x} should not be - * considered as an initial guess, as it is set to 0 in the initialization - * phase. - * @return A reference to {@code x} (shallow copy) updated with the - * solution. - * @throws NullArgumentException if one of the parameters is {@code null}. - * @throws NonSquareOperatorException if {@code a} or {@code m} is not - * square. - * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x} - * have dimensions inconsistent with {@code a}. - * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is - * {@code true}, and {@code a} or {@code m} is not self-adjoint. - * @throws IllConditionedOperatorException if {@code a} is ill-conditioned. - * @throws MaxCountExceededException at exhaustion of the iteration count, - * unless a custom {@link MaxCountExceededCallback callback} has been set at - * construction. - */ - @Override - public RealVector solveInPlace(final RealLinearOperator a, - final RealVector b, final RealVector x) - throws NullArgumentException, NonSquareOperatorException, - DimensionMismatchException, NonSelfAdjointOperatorException, - IllConditionedOperatorException, MaxCountExceededException { - return solveInPlace(a, null, b, x, false, 0.); - } - - /** - * Creates the event to be fired during the solution process. Unmodifiable - * views of the RHS vector, and the current estimate of the solution are - * returned by the created event. - * - * @param state Reference to the current state of this algorithm. - * @return The newly created event. - */ - private IterativeLinearSolverEvent createEvent(final State state) { - final RealVector bb = RealVector.unmodifiableRealVector(state.b); - - final IterativeLinearSolverEvent event; - event = new IterativeLinearSolverEvent(this) { - - @Override - public RealVector getRightHandSideVector() { - return bb; - } - - @Override - public RealVector getSolution() { - final int n = state.x.getDimension(); - final RealVector x = new ArrayRealVector(n); - state.refine(x); - return x; - } - }; - return event; - } -} +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.linear; + +import org.apache.commons.math.exception.DimensionMismatchException; +import org.apache.commons.math.exception.MaxCountExceededException; +import org.apache.commons.math.exception.NullArgumentException; +import org.apache.commons.math.exception.util.ExceptionContext; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.IterationManager; +import org.apache.commons.math.util.MathUtils; + +/** + *

+ * Implementation of the SYMMLQ iterative linear solver proposed by Paige and Saunders (1975). This implementation is + * largely based on the FORTRAN code by Pr. Michael A. Saunders, available here. + *

+ *

+ * SYMMLQ is designed to solve the system of linear equations A · x = b + * where A is an n × n self-adjoint linear operator (defined as a + * {@link RealLinearOperator}), and b is a given vector. The operator A is not + * required to be positive definite. If A is known to be definite, the method of + * conjugate gradients might be preferred, since it will require about the same + * number of iterations as SYMMLQ but slightly less work per iteration. + *

+ *

+ * SYMMLQ is designed to solve the system (A - shift · I) · x = b, + * where shift is a specified scalar value. If shift and b are suitably chosen, + * the computed vector x may approximate an (unnormalized) eigenvector of A, as + * in the methods of inverse iteration and/or Rayleigh-quotient iteration. + * Again, the linear operator (A - shift · I) need not be positive + * definite (but must be self-adjoint). The work per iteration is very + * slightly less if shift = 0. + *

+ *

Peconditioning

+ *

+ * Preconditioning may reduce the number of iterations required. The solver is + * provided with a positive definite preconditioner M = C · C^T + * that is known to approximate (A - shift · I) in some sense, while + * systems of the form M · y = x can be solved efficiently. Then SYMMLQ + * will implicitly solve the system of equations P · (A - shift · + * I) · P^T · xhat = P · b, i.e. Ahat · + * xhat = bhat, where P = C^-1, Ahat = P · (A - shift · + * I) · P^T, bhat = P · b, and return the solution x = + * P^T · xhat. The associated residual is rhat = bhat - Ahat + * · xhat = P · [b - (A - shift · I) · x] = P + * · r. + *

+ *

Default stopping criterion

+ *

+ * A default stopping criterion is implemented. The iterations stop when || rhat + * || ≤ δ || Ahat || || xhat ||, where xhat is the current estimate of + * the solution of the transformed system, rhat the current estimate of the + * corresponding residual, and δ a user-specified tolerance. + *

+ *

Iteration count

+ *

+ * In the present context, an iteration should be understood as one evaluation + * of the matrix-vector product A · x. The initialization phase therefore + * counts as one iteration. If the user requires checks on the symmetry of A, + * this entails one further matrix-vector product by iteration. This further + * product is not accounted for in the iteration count. In other words, + * the number of iterations required to reach convergence will be identical, + * whether checks have been required or not. + *

+ *

+ * The present definition of the iteration count differs from that adopted in + * the original FOTRAN code, where the initialization phase was not + * taken into account. + *

+ *

Exception context

+ *

+ * Besides standard {@link DimensionMismatchException}, this class might throw + * {@link NonSelfAdjointOperatorException} if the linear operator or the + * preconditioner are not symmetric. In this case, the {@link ExceptionContext} + * provides more information + *

key {@code "operator"} points to the offending linear operator, say L,
key {@code "vector1"} points to the first offending vector, say x, + *
key {@code "vector2"} points to the second offending vector, say y, such + * that x^T · L · y ≠ y^T · L + * · x (within a certain accuracy).

+ *

+ * {@link NonPositiveDefiniteOperatorException} might also be thrown in case the + * preconditioner is not positive definite. The relevant keys to the + * {@link ExceptionContext} are + *

key {@code "operator"}, which points to the offending linear operator, + * say L,
key {@code "vector"}, which points to the offending vector, say x, such + * that x^T · L · x < 0.

+ *

References

+ *

Paige and Saunders (1975): C. C. Paige and M. A. Saunders, + * Solution of Sparse Indefinite Systems of Linear Equations, SIAM + * Journal on Numerical Analysis 12(4): 617-629, 1975

+ * + * @version $Id$ + * @since 3.0 + */ +public class SymmLQ + extends PreconditionedIterativeLinearSolver { + + /* + * IMPLEMENTATION NOTES + * -------------------- + * The implementation follows as closely as possible the notations of Paige + * and Saunders (1975). Attention must be paid to the fact that some + * quantities which are relevant to iteration k can only be computed in + * iteration (k+1). Therefore, minute attention must be paid to the index of + * each state variable of this algorithm. + * + * 1. Preconditioning + * --------------- + * The Lanczos iterations associated with Ahat and bhat read + * beta[1] = |P . b| + * v[1] = P.b / beta[1] + * beta[k+1] * v[k+1] = Ahat * v[k] - alpha[k] * v[k] - beta[k] * v[k-1] + * = P * (A - shift * I) * P' * v[k] - alpha[k] * v[k] + * - beta[k] * v[k-1] + * Multiplying both sides by P', we get + * beta[k+1] * (P' * v)[k+1] = M^(-1) * (A - shift * I) * (P' * v)[k] + * - alpha[k] * (P' * v)[k] + * - beta[k] * (P' * v[k-1]), + * and + * alpha[k+1] = v[k+1]' * Ahat * v[k+1] + * = v[k+1]' * P * (A - shift * I) * P' * v[k+1] + * = (P' * v)[k+1]' * (A - shift * I) * (P' * v)[k+1]. + * + * In other words, the Lanczos iterations are unchanged, except for the fact + * that we really compute (P' * v) instead of v. It can easily be checked + * that all other formulas are unchanged. It must be noted that P is never + * explicitly used, only matrix-vector products involving M^(-1) are + * invoked. + * + * 2. Accounting for the shift parameter + * ---------------------------------- + * Is trivial: each time A.operate(x) is invoked, one must subtract shift * x + * to the result. + * + * 3. Accounting for the goodb flag + * ----------------------------- + * When goodb is set to true, the component of xL along b is computed + * separately. From Page and Saunders (1975), equation (5.9), we have + * wbar[k+1] = s[k] * wbar[k] - c[k] * v[k+1], + * wbar[1] = v[1]. + * Introducing wbar2[k] = wbar[k] - s[1] * ... * s[k-1] * v[1], it can + * easily be verified by induction that what follows the same recursive + * relation + * wbar2[k+1] = s[k] * wbar2[k] - c[k] * v[k+1], + * wbar2[1] = 0, + * and we then have + * w[k] = c[k] * wbar2[k] + s[k] * v[k+1] + * + s[1] * ... * s[k-1] * c[k] * v[1]. + * Introducing w2[k] = w[k] - s[1] * ... * s[k-1] * c[k] * v[1], we find, + * from (5.10) + * xL[k] = zeta[1] * w[1] + ... + zeta[k] * w[k] + * = zeta[1] * w2[1] + ... + zeta[k] * w2[k] + * + (s[1] * c[2] * zeta[2] + ... + * + s[1] * ... * s[k-1] * c[k] * zeta[k]) * v[1] + * = xL2[k] + bstep[k] * v[1], + * where xL2[k] is defined by + * xL2[0] = 0, + * xL2[k+1] = xL2[k] + zeta[k+1] * w2[k+1], + * and bstep is defined by + * bstep[1] = 0, + * bstep[k] = bstep[k-1] + s[1] * ... * s[k-1] * c[k] * zeta[k]. + * We also have, from (5.11) + * xC[k] = xL[k-1] + zbar[k] * wbar[k] + * = xL2[k-1] + zbar[k] * wbar2[k] + * + (bstep[k-1] + s[1] * ... * s[k-1] * zbar[k]) * v[1]. + */ + + /** + *

+ * A simple container holding the non-final variables used in the + * iterations. Making the current state of the solver visible from the + * outside is necessary, because during the iterations, {@code x} does not + * exactly hold the current estimate of the solution. Indeed, + * {@code x} needs in general to be moved from the LQ point to the CG point. + * Besides, additional upudates must be carried out in case {@code goodb} is + * set to {@code true}. + *

+ *

+ * In all subsequent comments, the description of the state variables refer + * to their value after a call to {@link #update()}. In these comments, k is + * the current number of evaluations of matrix-vector products. + *

+ */ + private class State { + + /** Reference to the linear operator. */ + private final RealLinearOperator a; + + /** Reference to the right-hand side vector. */ + private final RealVector b; + + /** The value of beta[k+1]. */ + private double beta; + + /** The value of beta[1]. */ + private double beta1; + + /** The value of bstep[k-1]. */ + private double bstep; + + /** The estimate of the norm of P * rC[k]. */ + private double cgnorm; + + /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */ + private double dbar; + + /** + * The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the + * initial code. + */ + private double gammaZeta; + + /** The value of gbar[k]. */ + private double gbar; + + /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */ + private double gmax; + + /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */ + private double gmin; + + /** Copy of the {@code goodb} parameter. */ + private final boolean goodb; + + /** {@code true} if the default convergence criterion is verified. */ + private boolean hasConverged; + + /** The estimate of the norm of P * rL[k-1]. */ + private double lqnorm; + + /** Reference to the preconditioner. */ + private final InvertibleRealLinearOperator m; + + /** + * The value of (-eps[k+1] * zeta[k-1]). Was called {@code rhs2} in the + * initial code. + */ + private double minusEpsZeta; + + /** The value of M^(-1) * b. */ + private final RealVector mSolveB; + + /** The value of beta[k]. */ + private double oldb; + + /** The value of beta[k] * M * P' * v[k]. */ + private RealVector r1; + + /** The value of beta[k+1] * M * P' * v[k+1]. */ + private RealVector r2; + + /** Copy of the {@code shift} parameter. */ + private final double shift; + + /** The value of s[1] * ... * s[k-1]. */ + private double snprod; + + /** + * An estimate of the square of the norm of A * V[k], based on Paige and + * Saunders (1975), equation (3.3). + */ + private double tnorm; + + /** + * The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] * + * v[1]) if {@code goodb} is {@code true}. Was called {@code w} in the + * initial code. + */ + private RealVector wbar; + + /** + * A reference to the vector to be updated with the solution. Contains + * the value of xL[k-1] if {@code goodb} is {@code false}, (xL[k-1] - + * bstep[k-1] * v[1]) otherwise. + */ + private final RealVector x; + + /** The value of beta[k+1] * P' * v[k+1]. */ + private RealVector y; + + /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */ + private double ynorm2; + + /** + * Creates and inits to k = 1 a new instance of this class. + * + * @param a Linear operator A of the system. + * @param m Preconditioner (can be {@code null}). + * @param b Right-hand side vector. + * @param x Vector to be updated with the solution. {@code x} should not + * be considered as an initial guess, as it is set to 0 in the + * initialization phase. + * @param goodb Usually {@code false}, except if {@code x} is expected + * to contain a large multiple of {@code b}. + * @param shift The amount to be subtracted to all diagonal elements of + * A. + */ + public State(final RealLinearOperator a, + final InvertibleRealLinearOperator m, final RealVector b, + final RealVector x, final boolean goodb, final double shift) { + this.a = a; + this.m = m; + this.b = b; + this.x = x; + this.goodb = goodb; + this.shift = shift; + this.mSolveB = m == null ? b : m.solve(b); + this.hasConverged = false; + init(); + } + + /** + * Move to the CG point if it seems better. In this version of SYMMLQ, + * the convergence tests involve only cgnorm, so we're unlikely to stop + * at an LQ point, except if the iteration limit interferes. + * + * @param xRefined Vector to be updated with the refined value of x. + */ + public void refine(final RealVector xRefined) { + final int n = this.x.getDimension(); + if (lqnorm < cgnorm) { + if (!goodb) { + xRefined.setSubVector(0, this.x); + } else { + final double step = bstep / beta1; + for (int i = 0; i < n; i++) { + final double bi = mSolveB.getEntry(i); + final double xi = this.x.getEntry(i); + xRefined.setEntry(i, xi + step * bi); + } + } + } else { + final double anorm = FastMath.sqrt(tnorm); + final double diag = gbar == 0. ? anorm * MACH_PREC : gbar; + final double zbar = gammaZeta / diag; + final double step = (bstep + snprod * zbar) / beta1; + // ynorm = FastMath.sqrt(ynorm2 + zbar * zbar); + if (!goodb) { + for (int i = 0; i < n; i++) { + final double xi = this.x.getEntry(i); + final double wi = wbar.getEntry(i); + xRefined.setEntry(i, xi + zbar * wi); + } + } else { + for (int i = 0; i < n; i++) { + final double xi = this.x.getEntry(i); + final double wi = wbar.getEntry(i); + final double bi = mSolveB.getEntry(i); + xRefined.setEntry(i, xi + zbar * wi + step * bi); + } + } + } + } + + /** + * Performs the initial phase of the SYMMLQ algorithm. On return, the + * value of the state variables of {@code this} object correspond to k = + * 1. + */ + private void init() { + this.x.set(0.); + /* + * Set up y for the first Lanczos vector. y and beta1 will be zero + * if b = 0. + */ + this.r1 = this.b.copy(); + this.y = this.m == null ? this.b.copy() : this.m.solve(this.r1); + if ((this.m != null) && check) { + checkSymmetry(this.m, this.r1, this.y, this.m.solve(this.y)); + } + + this.beta1 = this.r1.dotProduct(this.y); + if (this.beta1 < 0.) { + throwNPDLOException(this.m, this.y); + } + if (this.beta1 == 0.) { + /* If b = 0 exactly, stop with x = 0. */ + return; + } + this.beta1 = FastMath.sqrt(this.beta1); + /* At this point + * r1 = b, + * y = M^(-1) * b, + * beta1 = beta[1]. + */ + final RealVector v = this.y.mapMultiply(1. / this.beta1); + this.y = this.a.operate(v); + if (check) { + checkSymmetry(this.a, v, this.y, this.a.operate(this.y)); + } + /* + * Set up y for the second Lanczos vector. y and beta will be zero + * or very small if b is an eigenvector. + */ + daxpy(-this.shift, v, this.y); + final double alpha = v.dotProduct(this.y); + daxpy(-alpha / this.beta1, this.r1, this.y); + /* + * At this point + * alpha = alpha[1] + * y = beta[2] * M * P' * v[2] + */ + /* Make sure r2 will be orthogonal to the first v. */ + final double vty = v.dotProduct(this.y); + final double vtv = v.dotProduct(v); + daxpy(-vty / vtv, v, this.y); + this.r2 = this.y.copy(); + if (this.m != null) { + this.y = this.m.solve(this.r2); + } + this.oldb = this.beta1; + this.beta = this.r2.dotProduct(this.y); + if (this.beta < 0.) { + throwNPDLOException(this.m, this.y); + } + this.beta = FastMath.sqrt(this.beta); + /* + * At this point + * oldb = beta[1] + * beta = beta[2] + * y = beta[2] * P' * v[2] + * r2 = beta[2] * M * P' * v[2] + */ + this.cgnorm = this.beta1; + this.gbar = alpha; + this.dbar = this.beta; + this.gammaZeta = this.beta1; + this.minusEpsZeta = 0.; + this.bstep = 0.; + this.snprod = 1.; + this.tnorm = alpha * alpha + this.beta * this.beta; + this.ynorm2 = 0.; + this.gmax = FastMath.abs(alpha) + MACH_PREC; + this.gmin = this.gmax; + + if (this.goodb) { + this.wbar = new ArrayRealVector(this.a.getRowDimension()); + this.wbar.set(0.); + } else { + this.wbar = v; + } + updateNorms(); + } + + /** + * Performs the next iteration of the algorithm. The iteration count + * should be incremented prior to calling this method. On return, the + * value of the state variables of {@code this} object correspond to the + * current iteration count {@code k}. + */ + private void update() { + final RealVector v = y.mapMultiply(1. / beta); + y = a.operate(v); + daxpbypz(-shift, v, -beta / oldb, r1, y); + final double alpha = v.dotProduct(y); + /* + * At this point + * v = P' * v[k], + * y = (A - shift * I) * P' * v[k] - beta[k] * M * P' * v[k-1], + * alpha = v'[k] * P * (A - shift * I) * P' * v[k] + * - beta[k] * v[k]' * P * M * P' * v[k-1] + * = v'[k] * P * (A - shift * I) * P' * v[k] + * - beta[k] * v[k]' * v[k-1] + * = alpha[k]. + */ + daxpy(-alpha / beta, r2, y); + /* + * At this point + * y = (A - shift * I) * P' * v[k] - alpha[k] * M * P' * v[k] + * - beta[k] * M * P' * v[k-1] + * = M * P' * (P * (A - shift * I) * P' * v[k] -alpha[k] * v[k] + * - beta[k] * v[k-1]) + * = beta[k+1] * M * P' * v[k+1], + * from Paige and Saunders (1975), equation (3.2). + * + * WATCH-IT: the two following line work only because y is no longer + * updated up to the end of the present iteration, and is + * reinitialized at the beginning of the next iteration. + */ + r1 = r2; + r2 = y; + if (m != null) { + y = m.solve(r2); + } + oldb = beta; + beta = r2.dotProduct(y); + if (beta < 0.) { + throwNPDLOException(m, y); + } + beta = FastMath.sqrt(beta); + /* + * At this point + * r1 = beta[k] * M * P' * v[k], + * r2 = beta[k+1] * M * P' * v[k+1], + * y = beta[k+1] * P' * v[k+1], + * oldb = beta[k], + * beta = beta[k+1]. + */ + tnorm += alpha * alpha + oldb * oldb + beta * beta; + /* + * Compute the next plane rotation for Q. See Paige and Saunders + * (1975), equation (5.6), with + * gamma = gamma[k-1], + * c = c[k-1], + * s = s[k-1]. + */ + final double gamma = FastMath.sqrt(gbar * gbar + oldb * oldb); + final double c = gbar / gamma; + final double s = oldb / gamma; + /* + * The relations + * gbar[k] = s[k-1] * (-c[k-2] * beta[k]) - c[k-1] * alpha[k] + * = s[k-1] * dbar[k] - c[k-1] * alpha[k], + * delta[k] = c[k-1] * dbar[k] + s[k-1] * alpha[k], + * are not stated in Paige and Saunders (1975), but can be retrieved + * by expanding the (k, k-1) and (k, k) coefficients of the matrix in + * equation (5.5). + */ + final double deltak = c * dbar + s * alpha; + gbar = s * dbar - c * alpha; + final double eps = s * beta; + dbar = -c * beta; + final double zeta = gammaZeta / gamma; + /* + * At this point + * gbar = gbar[k] + * deltak = delta[k] + * eps = eps[k+1] + * dbar = dbar[k+1] + * zeta = zeta[k-1] + */ + final double zetaC = zeta * c; + final double zetaS = zeta * s; + final int n = x.getDimension(); + for (int i = 0; i < n; i++) { + final double xi = x.getEntry(i); + final double vi = v.getEntry(i); + final double wi = wbar.getEntry(i); + x.setEntry(i, xi + wi * zetaC + vi * zetaS); + wbar.setEntry(i, wi * s - vi * c); + } + /* + * At this point + * x = xL[k-1], + * ptwbar = P' wbar[k], + * see Paige and Saunders (1975), equations (5.9) and (5.10). + */ + bstep += snprod * c * zeta; + snprod *= s; + gmax = FastMath.max(gmax, gamma); + gmin = FastMath.min(gmin, gamma); + ynorm2 += zeta * zeta; + gammaZeta = minusEpsZeta - deltak * zeta; + minusEpsZeta = -eps * zeta; + /* + * At this point + * snprod = s[1] * ... * s[k-1], + * gmax = max(|alpha[1]|, gamma[1], ..., gamma[k-1]), + * gmin = min(|alpha[1]|, gamma[1], ..., gamma[k-1]), + * ynorm2 = zeta[1]^2 + ... + zeta[k-1]^2, + * gammaZeta = gamma[k] * zeta[k], + * minusEpsZeta = -eps[k+1] * zeta[k-1]. + * The relation for gammaZeta can be retrieved from Paige and + * Saunders (1975), equation (5.4a), last line of the vector + * gbar[k] * zbar[k] = -eps[k] * zeta[k-2] - delta[k] * zeta[k-1]. + */ + updateNorms(); + } + + /** + * Computes the norms of the residuals, and checks for convergence. + * Updates {@link #lqnorm} and {@link #cgnorm}. + */ + private void updateNorms() { + final double anorm = FastMath.sqrt(tnorm); + final double ynorm = FastMath.sqrt(ynorm2); + final double epsa = anorm * MACH_PREC; + final double epsx = anorm * ynorm * MACH_PREC; + final double epsr = anorm * ynorm * delta; + final double diag = gbar == 0. ? epsa : gbar; + lqnorm = FastMath.sqrt(gammaZeta * gammaZeta + + minusEpsZeta * minusEpsZeta); + final double qrnorm = snprod * beta1; + cgnorm = qrnorm * beta / FastMath.abs(diag); + + /* + * Estimate cond(A). In this version we look at the diagonals of L + * in the factorization of the tridiagonal matrix, T = L * Q. + * Sometimes, T[k] can be misleadingly ill-conditioned when T[k+1] + * is not, so we must be careful not to overestimate acond. + */ + final double acond; + if (lqnorm <= cgnorm) { + acond = gmax / gmin; + } else { + acond = gmax / FastMath.min(gmin, FastMath.abs(diag)); + } + if (acond * MACH_PREC >= 0.1) { + throw new IllConditionedOperatorException(acond); + } + if (beta1 <= epsx) { + /* + * x has converged to an eigenvector of A corresponding to the + * eigenvalue shift. + */ + throw new SingularOperatorException(); + } + hasConverged = (cgnorm <= epsx) || (cgnorm <= epsr); + } + } + + /** The cubic root of {@link #MACH_PREC}. */ + private static final double CBRT_MACH_PREC; + + /** The machine precision. */ + private static final double MACH_PREC; + + /** Key for the exception context. */ + private static final String OPERATOR = "operator"; + + /** Key for the exception context. */ + private static final String THRESHOLD = "threshold"; + + /** Key for the exception context. */ + private static final String VECTOR = "vector"; + + /** Key for the exception context. */ + private static final String VECTOR1 = "vector1"; + + /** Key for the exception context. */ + private static final String VECTOR2 = "vector2"; + + /** {@code true} if symmetry of matrix and conditioner must be checked. */ + private final boolean check; + + /** + * The value of the custom tolerance δ for the default stopping + * criterion. + */ + private final double delta; + + /** + * Creates a new instance of this class, with default + * stopping criterion. + * + * @param maxIterations Maximum number of iterations. + * @param delta δ parameter for the default stopping criterion. + * @param check {@code true} if self-adjointedness of both matrix and + * preconditioner should be checked. This entails an extra matrix-vector + * product at each iteration. + */ + public SymmLQ(final int maxIterations, final double delta, + final boolean check) { + super(maxIterations); + this.delta = delta; + this.check = check; + } + + /** + * Creates a new instance of this class, with default + * stopping criterion and custom iteration manager. + * + * @param manager Custom iteration manager. + * @param delta δ parameter for the default stopping criterion. + * @param check {@code true} if self-adjointedness of both matrix and + * preconditioner should be checked. This entails an extra matrix-vector + * product at each iteration. + */ + public SymmLQ(final IterationManager manager, final double delta, + final boolean check) { + super(manager); + this.delta = delta; + this.check = check; + } + + static { + MACH_PREC = Math.ulp(1.); + CBRT_MACH_PREC = Math.cbrt(MACH_PREC); + } + + /** + * Performs a symmetry check on the specified linear operator, and throws an + * exception in case this check fails. Given a linear operator L, and a + * vector x, this method checks that x' L y = y' L x (within a given + * accuracy), where y = L x. + * + * @param l The linear operator L. + * @param x The candidate vector x. + * @param y The candidate vector y = L x. + * @param z The vector z = L y. + * @throws NonSelfAdjointOperatorException when the test fails. + */ + private static void checkSymmetry(final RealLinearOperator l, + final RealVector x, final RealVector y, + final RealVector z) + throws NonSelfAdjointOperatorException { + final double s = y.dotProduct(y); + final double t = x.dotProduct(z); + final double epsa = (s + MACH_PREC) * CBRT_MACH_PREC; + if (FastMath.abs(s - t) > epsa) { + final NonSelfAdjointOperatorException e; + e = new NonSelfAdjointOperatorException(); + final ExceptionContext context = e.getContext(); + context.setValue(OPERATOR, l); + context.setValue(VECTOR1, x); + context.setValue(VECTOR2, y); + context.setValue(THRESHOLD, Double.valueOf(epsa)); + throw e; + } + } + + /** + * A BLAS-like function, for the operation z ← a · x + b + * · y + z. This is for internal use only: no dimension checks are + * provided. + * + * @param a The scalar by which {@code x} is to be multiplied. + * @param x The first vector to be added to {@code z}. + * @param b The scalar by which {@code y} is to be multiplied. + * @param y The second vector to be added to {@code z}. + * @param z The vector to be incremented. + */ + private static void daxpbypz(final double a, final RealVector x, + final double b, final RealVector y, + final RealVector z) { + final int n = z.getDimension(); + for (int i = 0; i < n; i++) { + final double zi; + zi = a * x.getEntry(i) + b * y.getEntry(i) + z.getEntry(i); + z.setEntry(i, zi); + } + } + + /** + * A clone of the BLAS {@code DAXPY} function, which carries out the + * operation y ← a · x + y. This is for internal use only: no + * dimension checks are provided. + * + * @param a The scalar by which {@code x} is to be multiplied. + * @param x The vector to be added to {@code y}. + * @param y The vector to be incremented. + */ + private static void daxpy(final double a, final RealVector x, + final RealVector y) { + final int n = x.getDimension(); + for (int i = 0; i < n; i++) { + y.setEntry(i, a * x.getEntry(i) + y.getEntry(i)); + } + } + + /** + * Throws a new {@link NonPositiveDefiniteOperatorException} with + * appropriate context. + * + * @param l The offending linear operator. + * @param v The offending vector. + * @throws NonPositiveDefiniteOperatorException in any circumstances. + */ + private static void throwNPDLOException(final RealLinearOperator l, + final RealVector v) + throws NonPositiveDefiniteOperatorException { + final NonPositiveDefiniteOperatorException e; + e = new NonPositiveDefiniteOperatorException(); + final ExceptionContext context = e.getContext(); + context.setValue(OPERATOR, l); + context.setValue(VECTOR, v); + throw e; + } + + /** + * Returns {@code true} if symmetry of the matrix, and symmetry as well as + * positive definiteness of the preconditioner should be checked. + * + * @return {@code true} if the tests are to be performed. + */ + public final boolean getCheck() { + return check; + } + + /** + * Returns an estimate of the solution to the linear system A · x = + * b. + * [... 412 lines stripped ...]