From 17f788387d8ff046d55c5b0f2eb72c7a98b54ff3 Mon Sep 17 00:00:00 2001 From: Sebastien Brisard Date: Sat, 22 Oct 2011 06:36:31 +0000 Subject: [PATCH] Implementation of the SYMMLQ iterative linear solver, based on Pr. Saunders FORTRAN impl. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1187657 13f79535-47bb-0310-9956-ffa450edef68 --- .../apache/commons/math/linear/SymmLQ.java | 1223 +++++++++++++++++ .../commons/math/linear/SymmLQTest.java | 641 +++++++++ 2 files changed, 1864 insertions(+) create mode 100644 src/main/java/org/apache/commons/math/linear/SymmLQ.java create mode 100644 src/test/java/org/apache/commons/math/linear/SymmLQTest.java diff --git a/src/main/java/org/apache/commons/math/linear/SymmLQ.java b/src/main/java/org/apache/commons/math/linear/SymmLQ.java new file mode 100644 index 000000000..955763073 --- /dev/null +++ b/src/main/java/org/apache/commons/math/linear/SymmLQ.java @@ -0,0 +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 · CT + * 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) · PT · xhat = P · b, i.e. Ahat · + * xhat = bhat, where P = C-1, Ahat = P · (A - shift · + * I) · PT, bhat = P · b, and return the solution x = + * PT · 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 + *

+ *

+ *

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

+ *

+ *

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 bT · A · b / + * (bT · 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. + *

+ * 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}. + *

+ *

+ * {@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 bT · A · b / + * (bT · 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 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. + *

+ * 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 bT · A · b / + * (bT · 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 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; + } +} diff --git a/src/test/java/org/apache/commons/math/linear/SymmLQTest.java b/src/test/java/org/apache/commons/math/linear/SymmLQTest.java new file mode 100644 index 000000000..92d0fc148 --- /dev/null +++ b/src/test/java/org/apache/commons/math/linear/SymmLQTest.java @@ -0,0 +1,641 @@ +/* + * 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.util.FastMath; +import org.apache.commons.math.util.IterationEvent; +import org.apache.commons.math.util.IterationListener; +import org.junit.Assert; +import org.junit.Test; + +public class SymmLQTest { + + public void saundersTest(final int n, final boolean goodb, + final boolean precon, final double shift, + final double pertbn) { + final RealLinearOperator a = new RealLinearOperator() { + + @Override + public RealVector operate(final RealVector x) { + if (x.getDimension() != n) { + throw new DimensionMismatchException(x.getDimension(), n); + } + final double[] y = new double[n]; + for (int i = 0; i < n; i++) { + y[i] = (i + 1) * 1.1 / n * x.getEntry(i); + } + return new ArrayRealVector(y, false); + } + + @Override + public int getRowDimension() { + return n; + } + + @Override + public int getColumnDimension() { + return n; + } + }; + final double shiftm = shift; + final double pertm = FastMath.abs(pertbn); + final InvertibleRealLinearOperator m; + if (precon) { + m = new InvertibleRealLinearOperator() { + + @Override + public RealVector operate(final RealVector x) { + if (x.getDimension() != n) { + throw new DimensionMismatchException(x.getDimension(), + n); + } + final double[] y = new double[n]; + for (int i = 0; i < n; i++) { + double d = (i + 1) * 1.1 / n; + d = FastMath.abs(d - shiftm); + if (i % 10 == 0) { + d += pertm; + } + y[i] = d * x.getEntry(i); + } + return new ArrayRealVector(y, false); + } + + @Override + public int getRowDimension() { + return n; + } + + @Override + public int getColumnDimension() { + return n; + } + + @Override + public RealVector solve(final RealVector b) { + if (b.getDimension() != n) { + throw new DimensionMismatchException(b.getDimension(), + n); + } + final double[] x = new double[n]; + for (int i = 0; i < n; i++) { + double d = (i + 1) * 1.1 / n; + d = FastMath.abs(d - shiftm); + if (i % 10 == 0) { + d += pertm; + } + x[i] = b.getEntry(i) / d; + } + return new ArrayRealVector(x, false); + } + }; + } else { + m = null; + } + final RealVector xtrue = new ArrayRealVector(n); + for (int i = 0; i < n; i++) { + xtrue.setEntry(i, n - i); + } + final RealVector b = a.operate(xtrue); + b.combineToSelf(1.0, -shift, xtrue); + final SymmLQ solver = new SymmLQ(2 * n, 1E-12, true); + final RealVector x = solver.solve(a, m, b, goodb, shift); + final RealVector y = a.operate(x); + final RealVector r1 = new ArrayRealVector(n); + for (int i = 0; i < n; i++) { + final double bi = b.getEntry(i); + final double yi = y.getEntry(i); + final double xi = x.getEntry(i); + r1.setEntry(i, bi - yi + shift * xi); + } + final double enorm = x.subtract(xtrue).getNorm() / xtrue.getNorm(); + final double etol = 1E-5; + Assert.assertTrue("enorm=" + + enorm + + ", " + + solver.getIterationManager() + .getIterations(), + enorm <= etol); + } + + @Test + public void testSolveSaunders1() { + saundersTest(1, false, false, 0., 0.); + } + + @Test + public void testSolveSaunders2() { + saundersTest(2, false, false, 0., 0.); + } + + @Test + public void testSolveSaunders3() { + saundersTest(1, false, true, 0., 0.); + } + + @Test + public void testSolveSaunders4() { + saundersTest(2, false, true, 0., 0.); + } + + @Test + public void testSolveSaunders5() { + saundersTest(5, false, true, 0., 0.); + } + + @Test + public void testSolveSaunders6() { + saundersTest(5, false, true, 0.25, 0.); + } + + @Test + public void testSolveSaunders7() { + saundersTest(50, false, false, 0., 0.); + } + + @Test + public void testSolveSaunders8() { + saundersTest(50, false, false, 0.25, 0.); + } + + @Test + public void testSolveSaunders9() { + saundersTest(50, false, true, 0., 0.10); + } + + @Test + public void testSolveSaunders10() { + saundersTest(50, false, true, 0.25, 0.10); + } + + @Test + public void testSolveSaunders11() { + saundersTest(1, true, false, 0., 0.); + } + + @Test + public void testSolveSaunders12() { + saundersTest(2, true, false, 0., 0.); + } + + @Test + public void testSolveSaunders13() { + saundersTest(1, true, true, 0., 0.); + } + + @Test + public void testSolveSaunders14() { + saundersTest(2, true, true, 0., 0.); + } + + @Test + public void testSolveSaunders15() { + saundersTest(5, true, true, 0., 0.); + } + + @Test + public void testSolveSaunders16() { + saundersTest(5, true, true, 0.25, 0.); + } + + @Test + public void testSolveSaunders17() { + saundersTest(50, true, false, 0., 0.); + } + + @Test + public void testSolveSaunders18() { + saundersTest(50, true, false, 0.25, 0.); + } + + @Test + public void testSolveSaunders19() { + saundersTest(50, true, true, 0., 0.10); + } + + @Test + public void testSolveSaunders20() { + saundersTest(50, true, true, 0.25, 0.10); + } + + @Test(expected = NonSquareOperatorException.class) + public void testNonSquareOperator() { + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 3); + final IterativeLinearSolver solver; + solver = new SymmLQ(10, 0., false); + final ArrayRealVector b = new ArrayRealVector(a.getRowDimension()); + final ArrayRealVector x = new ArrayRealVector(a.getColumnDimension()); + solver.solve(a, b, x); + } + + @Test(expected = DimensionMismatchException.class) + public void testDimensionMismatchRightHandSide() { + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(3, 3); + final IterativeLinearSolver solver; + solver = new SymmLQ(10, 0., false); + final ArrayRealVector b = new ArrayRealVector(2); + solver.solve(a, b); + } + + @Test(expected = DimensionMismatchException.class) + public void testDimensionMismatchSolution() { + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(3, 3); + final IterativeLinearSolver solver; + solver = new SymmLQ(10, 0., false); + final ArrayRealVector b = new ArrayRealVector(3); + final ArrayRealVector x = new ArrayRealVector(2); + solver.solve(a, b, x); + } + + @Test + public void testUnpreconditionedSolution() { + final int n = 5; + final int maxIterations = 100; + final RealLinearOperator a = new HilbertMatrix(n); + final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n); + final IterativeLinearSolver solver; + solver = new SymmLQ(maxIterations, 1E-10, true); + final RealVector b = new ArrayRealVector(n); + for (int j = 0; j < n; j++) { + b.set(0.); + b.setEntry(j, 1.); + final RealVector x = solver.solve(a, b); + for (int i = 0; i < n; i++) { + final double actual = x.getEntry(i); + final double expected = ainv.getEntry(i, j); + final double delta = 1E-6 * Math.abs(expected); + final String msg = String.format("entry[%d][%d]", i, j); + Assert.assertEquals(msg, expected, actual, delta); + } + } + } + + @Test + public void testUnpreconditionedInPlaceSolutionWithInitialGuess() { + final int n = 5; + final int maxIterations = 100; + final RealLinearOperator a = new HilbertMatrix(n); + final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n); + final IterativeLinearSolver solver; + solver = new SymmLQ(maxIterations, 1E-10, true); + final RealVector b = new ArrayRealVector(n); + for (int j = 0; j < n; j++) { + b.set(0.); + b.setEntry(j, 1.); + final RealVector x0 = new ArrayRealVector(n); + x0.set(1.); + final RealVector x = solver.solveInPlace(a, b, x0); + Assert.assertSame("x should be a reference to x0", x0, x); + for (int i = 0; i < n; i++) { + final double actual = x.getEntry(i); + final double expected = ainv.getEntry(i, j); + final double delta = 1E-6 * Math.abs(expected); + final String msg = String.format("entry[%d][%d)", i, j); + Assert.assertEquals(msg, expected, actual, delta); + } + } + } + + @Test + public void testUnpreconditionedSolutionWithInitialGuess() { + final int n = 5; + final int maxIterations = 100; + final RealLinearOperator a = new HilbertMatrix(n); + final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n); + final IterativeLinearSolver solver; + solver = new SymmLQ(maxIterations, 1E-10, true); + final RealVector b = new ArrayRealVector(n); + for (int j = 0; j < n; j++) { + b.set(0.); + b.setEntry(j, 1.); + final RealVector x0 = new ArrayRealVector(n); + x0.set(1.); + final RealVector x = solver.solve(a, b, x0); + Assert.assertNotSame("x should not be a reference to x0", x0, x); + for (int i = 0; i < n; i++) { + final double actual = x.getEntry(i); + final double expected = ainv.getEntry(i, j); + final double delta = 1E-6 * Math.abs(expected); + final String msg = String.format("entry[%d][%d]", i, j); + Assert.assertEquals(msg, expected, actual, delta); + Assert.assertEquals(msg, x0.getEntry(i), 1., Math.ulp(1.)); + } + } + } + + @Test(expected = NonSquareOperatorException.class) + public void testNonSquarePreconditioner() { + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2); + final InvertibleRealLinearOperator m; + m = new InvertibleRealLinearOperator() { + + @Override + public RealVector operate(final RealVector x) { + throw new UnsupportedOperationException(); + } + + @Override + public int getRowDimension() { + return 2; + } + + @Override + public int getColumnDimension() { + return 3; + } + + @Override + public RealVector solve(final RealVector b) { + throw new UnsupportedOperationException(); + } + }; + final PreconditionedIterativeLinearSolver solver; + solver = new SymmLQ(10, 0., false); + final ArrayRealVector b = new ArrayRealVector(a.getRowDimension()); + solver.solve(a, m, b); + } + + @Test(expected = DimensionMismatchException.class) + public void testMismatchedOperatorDimensions() { + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2); + final InvertibleRealLinearOperator m; + m = new InvertibleRealLinearOperator() { + + @Override + public RealVector operate(final RealVector x) { + throw new UnsupportedOperationException(); + } + + @Override + public int getRowDimension() { + return 3; + } + + @Override + public int getColumnDimension() { + return 3; + } + + @Override + public RealVector solve(final RealVector b) { + throw new UnsupportedOperationException(); + } + }; + final PreconditionedIterativeLinearSolver solver; + solver = new SymmLQ(10, 0d, false); + final ArrayRealVector b = new ArrayRealVector(a.getRowDimension()); + solver.solve(a, m, b); + } + + @Test(expected = NonPositiveDefiniteOperatorException.class) + public void testNonPositiveDefinitePreconditioner() { + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2); + a.setEntry(0, 0, 1d); + a.setEntry(0, 1, 2d); + a.setEntry(1, 0, 3d); + a.setEntry(1, 1, 4d); + final InvertibleRealLinearOperator m; + m = new InvertibleRealLinearOperator() { + + @Override + public RealVector operate(final RealVector x) { + final ArrayRealVector y = new ArrayRealVector(2); + y.setEntry(0, -x.getEntry(0)); + y.setEntry(1, -x.getEntry(1)); + return y; + } + + @Override + public int getRowDimension() { + return 2; + } + + @Override + public int getColumnDimension() { + return 2; + } + + @Override + public RealVector solve(final RealVector b) { + final ArrayRealVector x = new ArrayRealVector(2); + x.setEntry(0, -b.getEntry(0)); + x.setEntry(1, -b.getEntry(1)); + return x; + } + }; + final PreconditionedIterativeLinearSolver solver; + solver = new SymmLQ(10, 0d, true); + final ArrayRealVector b = new ArrayRealVector(2); + b.setEntry(0, -1d); + b.setEntry(1, -1d); + solver.solve(a, m, b); + } + + @Test + public void testPreconditionedSolution() { + final int n = 8; + final int maxIterations = 100; + final RealLinearOperator a = new HilbertMatrix(n); + final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n); + final InvertibleRealLinearOperator m = JacobiPreconditioner.create(a); + final PreconditionedIterativeLinearSolver solver; + solver = new SymmLQ(maxIterations, 1E-15, true); + final RealVector b = new ArrayRealVector(n); + for (int j = 0; j < n; j++) { + b.set(0.); + b.setEntry(j, 1.); + final RealVector x = solver.solve(a, m, b); + for (int i = 0; i < n; i++) { + final double actual = x.getEntry(i); + final double expected = ainv.getEntry(i, j); + final double delta = 1E-6 * Math.abs(expected); + final String msg = String.format("coefficient (%d, %d)", i, j); + Assert.assertEquals(msg, expected, actual, delta); + } + } + } + + @Test + public void testPreconditionedSolution2() { + final int n = 100; + final int maxIterations = 100000; + final Array2DRowRealMatrix a = new Array2DRowRealMatrix(n, n); + double daux = 1.; + for (int i = 0; i < n; i++) { + a.setEntry(i, i, daux); + daux *= 1.2; + for (int j = i + 1; j < n; j++) { + if (i == j) { + } else { + final double value = 1.0; + a.setEntry(i, j, value); + a.setEntry(j, i, value); + } + } + } + final InvertibleRealLinearOperator m = JacobiPreconditioner.create(a); + final PreconditionedIterativeLinearSolver prec; + final IterativeLinearSolver unprec; + prec = new SymmLQ(maxIterations, 1E-15, true); + unprec = new SymmLQ(maxIterations, 1E-15, true); + final RealVector b = new ArrayRealVector(n); + final String pattern = "preconditioned SymmLQ (%d iterations) should" + + " have been faster than unpreconditioned (%d iterations)"; + String msg; + for (int j = 0; j < 1; j++) { + b.set(0.); + b.setEntry(j, 1.); + final RealVector px = prec.solve(a, m, b); + final RealVector x = unprec.solve(a, b); + final int npcg = prec.getIterationManager().getIterations(); + final int ncg = unprec.getIterationManager().getIterations(); + msg = String.format(pattern, npcg, ncg); + Assert.assertTrue(msg, npcg < ncg); + for (int i = 0; i < n; i++) { + msg = String.format("row %d, column %d", i, j); + final double expected = x.getEntry(i); + final double actual = px.getEntry(i); + final double delta = 5E-5 * Math.abs(expected); + Assert.assertEquals(msg, expected, actual, delta); + } + } + } + + @Test + public void testEventManagement() { + final int n = 5; + final int maxIterations = 100; + final RealLinearOperator a = new HilbertMatrix(n); + final IterativeLinearSolver solver; + final int[] count = new int[] { + 0, 0, 0, 0 + }; + final IterationListener listener = new IterationListener() { + + public void initializationPerformed(final IterationEvent e) { + count[0] = 1; + count[1] = 0; + count[2] = 0; + count[3] = 0; + + } + + public void iterationPerformed(final IterationEvent e) { + ++count[2]; + } + + public void iterationStarted(final IterationEvent e) { + ++count[1]; + + } + + public void terminationPerformed(final IterationEvent e) { + ++count[3]; + } + }; + solver = new SymmLQ(maxIterations, 1E-10, true); + solver.getIterationManager().addIterationListener(listener); + final RealVector b = new ArrayRealVector(n); + for (int j = 0; j < n; j++) { + b.set(0.); + b.setEntry(j, 1.); + solver.solve(a, b); + String msg = String.format("column %d (initialization)", j); + Assert.assertEquals(msg, 1, count[0]); + msg = String.format("column %d (iterations started)", j); + Assert.assertEquals(msg, solver.getIterationManager() + .getIterations() - 1, count[1]); + msg = String.format("column %d (iterations performed)", j); + Assert.assertEquals(msg, solver.getIterationManager() + .getIterations() - 1, count[2]); + msg = String.format("column %d (finalization)", j); + Assert.assertEquals(msg, 1, count[3]); + } + } + + @Test(expected = NonSelfAdjointOperatorException.class) + public void testNonSelfAdjointOperator() { + final RealLinearOperator a; + a = new Array2DRowRealMatrix(new double[][] { + { + 1., 2., 3. + }, { + 2., 4., 5. + }, { + 2.999, 5., 6. + } + }); + final RealVector b; + b = new ArrayRealVector(new double[] { + 1., 1., 1. + }); + new SymmLQ(100, 1., true).solve(a, b); + } + + @Test(expected = NonSelfAdjointOperatorException.class) + public void testNonSelfAdjointPreconditioner() { + final RealLinearOperator a = new Array2DRowRealMatrix(new double[][] { + { + 1., 2., 3. + }, { + 2., 4., 5. + }, { + 3., 5., 6. + } + }); + final Array2DRowRealMatrix mMat; + mMat = new Array2DRowRealMatrix(new double[][] { + { + 1., 0., 1. + }, { + 0., 1., 0. + }, { + 0., 0., 1. + } + }); + final DecompositionSolver mSolver; + mSolver = new LUDecomposition(mMat).getSolver(); + final InvertibleRealLinearOperator m; + m = new InvertibleRealLinearOperator() { + + @Override + public RealVector operate(final RealVector x) { + return mMat.operate(x); + } + + @Override + public int getRowDimension() { + return mMat.getRowDimension(); + } + + @Override + public int getColumnDimension() { + return mMat.getColumnDimension(); + } + + @Override + public RealVector solve(final RealVector b) { + return mSolver.solve(b); + } + }; + final RealVector b = new ArrayRealVector(new double[] { + 1., 1., 1. + }); + new SymmLQ(100, 1., true).solve(a, m, b); + } +}