commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From celes...@apache.org
Subject svn commit: r1187657 - 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 06:36:31 GMT
Author: celestin
Date: Sat Oct 22 06:36:31 2011
New Revision: 1187657

URL: http://svn.apache.org/viewvc?rev=1187657&view=rev
Log:
Implementation of the SYMMLQ iterative linear solver, based on Pr. Saunders FORTRAN impl.

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java

Added: 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=1187657&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java Sat Oct 22 06:36:31 2011
@@ -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;
+
+/**
+ * <p>
+ * Implementation of the SYMMLQ iterative linear solver proposed by <a
+ * href="#PAIG1975">Paige and Saunders (1975)</a>. This implementation is
+ * largely based on the FORTRAN code by Pr. Michael A. Saunders, available <a
+ * href="http://www.stanford.edu/group/SOL/software/symmlq/f77/">here</a>.
+ * </p>
+ * <p>
+ * SYMMLQ is designed to solve the system of linear equations A &middot; x = b
+ * where A is an n &times; 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.
+ * </p>
+ * <p>
+ * SYMMLQ is designed to solve the system (A - shift &middot; I) &middot; 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 &middot; I) need not be positive
+ * definite (but <em>must</em> be self-adjoint). The work per iteration is very
+ * slightly less if shift = 0.
+ * </p>
+ * <h3>Peconditioning</h3>
+ * <p>
+ * Preconditioning may reduce the number of iterations required. The solver is
+ * provided with a positive definite preconditioner M = C &middot; C<sup>T</sup>
+ * that is known to approximate (A - shift &middot; I) in some sense, while
+ * systems of the form M &middot; y = x can be solved efficiently. Then SYMMLQ
+ * will implicitly solve the system of equations P &middot; (A - shift &middot;
+ * I) &middot; P<sup>T</sup> &middot; xhat = P &middot; b, i.e. Ahat &middot;
+ * xhat = bhat, where P = C<sup>-1</sup>, Ahat = P &middot; (A - shift &middot;
+ * I) &middot; P<sup>T</sup>, bhat = P &middot; b, and return the solution x =
+ * P<sup>T</sup> &middot; xhat. The associated residual is rhat = bhat - Ahat
+ * &middot; xhat = P &middot; [b - (A - shift &middot; I) &middot; x] = P
+ * &middot; r.
+ * </p>
+ * <h3><a id="stopcrit">Default stopping criterion</a></h3>
+ * <p>
+ * A default stopping criterion is implemented. The iterations stop when || rhat
+ * || &le; &delta; || Ahat || || xhat ||, where xhat is the current estimate of
+ * the solution of the transformed system, rhat the current estimate of the
+ * corresponding residual, and &delta; a user-specified tolerance.
+ * </p>
+ * <h3>Iteration count</h3>
+ * <p>
+ * In the present context, an iteration should be understood as one evaluation
+ * of the matrix-vector product A &middot; 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 <em>not</em> 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.
+ * </p>
+ * <p>
+ * The present definition of the iteration count differs from that adopted in
+ * the original FOTRAN code, where the initialization phase was <em>not</em>
+ * taken into account.
+ * </p>
+ * <h3><a id="context">Exception context</a></h3>
+ * <p>
+ * 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
+ * <ul>
+ * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
+ * <li>key {@code "vector1"} points to the first offending vector, say x,
+ * <li>key {@code "vector2"} points to the second offending vector, say y, such
+ * that x<sup>T</sup> &middot; L &middot; y &ne; y<sup>T</sup> &middot; L
+ * &middot; x (within a certain accuracy).</li>
+ * </ul>
+ * </p>
+ * <p>
+ * {@link NonPositiveDefiniteOperatorException} might also be thrown in case the
+ * preconditioner is not positive definite. The relevant keys to the
+ * {@link ExceptionContext} are
+ * <ul>
+ * <li>key {@code "operator"}, which points to the offending linear operator,
+ * say L,</li>
+ * <li>key {@code "vector"}, which points to the offending vector, say x, such
+ * that x<sup>T</sup> &middot; L &middot; x < 0.</li>
+ * </ul>
+ * </p>
+ * <h3>References</h3>
+ * <dl>
+ * <dt><a id="PAIG1975">Paige and Saunders (1975)</a></dt>
+ * <dd>C. C. Paige and M. A. Saunders, <a
+ * href="http://www.stanford.edu/group/SOL/software/symmlq/PS75.pdf"><em>
+ * Solution of Sparse Indefinite Systems of Linear Equations</em></a>, SIAM
+ * Journal on Numerical Analysis 12(4): 617-629, 1975</dd>
+ * </dl>
+ *
+ * @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].
+     */
+
+    /**
+     * <p>
+     * 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
+     * <em>exactly</em> 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}.
+     * </p>
+     * <p>
+     * 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.
+     * </p>
+     */
+    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 &delta; for the default stopping
+     * criterion.
+     */
+    private final double delta;
+
+    /**
+     * Creates a new instance of this class, with <a href="#stopcrit">default
+     * stopping criterion</a>.
+     *
+     * @param maxIterations Maximum number of iterations.
+     * @param delta &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 <a href="#stopcrit">default
+     * stopping criterion</a> and custom iteration manager.
+     *
+     * @param manager Custom iteration manager.
+     * @param delta &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 &larr; a &middot; x + b
+     * &middot; 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 &larr; a &middot; 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 &middot; 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
+     * &middot; I) &middot; x = b.
+     * <p>
+     * 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.
+     * </p>
+     * <p>
+     * {@code shift} should be zero if the system A &middot; x = b is to be
+     * solved. Otherwise, it could be an approximation to an eigenvalue of A,
+     * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
+     * (b<sup>T</sup> &middot; 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.
+     * </p>
+     *
+     * @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 &middot; 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 &middot; 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 &middot; I) &middot; x = b.
+     * <p>
+     * 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}.
+     * </p>
+     * <p>
+     * {@code shift} should be zero if the system A &middot; x = b is to be
+     * solved. Otherwise, it could be an approximation to an eigenvalue of A,
+     * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
+     * (b<sup>T</sup> &middot; 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.
+     * </p>
+     *
+     * @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 &middot; 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 &middot; 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
+     * &middot; I) &middot; x = b. The solution is computed in-place.
+     * <p>
+     * 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.
+     * </p>
+     * <p>
+     * {@code shift} should be zero if the system A &middot; x = b is to be
+     * solved. Otherwise, it could be an approximation to an eigenvalue of A,
+     * such as the Rayleigh quotient b<sup>T</sup> &middot; A &middot; b /
+     * (b<sup>T</sup> &middot; 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.
+     * </p>
+     *
+     * @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 &middot; 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;
+    }
+}

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java?rev=1187657&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java Sat Oct 22 06:36:31 2011
@@ -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);
+    }
+}



Mime
View raw message