Class SymmLQ.State

  • Enclosing class:
    SymmLQ

    private static class SymmLQ.State
    extends java.lang.Object

    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, x does not exactly hold the current estimate of the solution. Indeed, x needs in general to be moved from the LQ point to the CG point. Besides, additional upudates must be carried out in case goodb is set to true.

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

    • Field Summary

      Fields 
      Modifier and Type Field Description
      private RealLinearOperator a
      Reference to the linear operator.
      private RealVector b
      Reference to the right-hand side vector.
      private double beta
      The value of beta[k+1].
      private double beta1
      The value of beta[1].
      private boolean bIsNull
      The value of b == 0 (exact floating-point equality).
      private double bstep
      The value of bstep[k-1].
      (package private) static double CBRT_MACH_PREC
      The cubic root of MACH_PREC.
      private double cgnorm
      The estimate of the norm of P * rC[k].
      private boolean check
      true if symmetry of matrix and conditioner must be checked.
      private double dbar
      The value of dbar[k+1] = -beta[k+1] * c[k-1].
      private double delta
      The value of the custom tolerance δ for the default stopping criterion.
      private double gammaZeta
      The value of gamma[k] * zeta[k].
      private double gbar
      The value of gbar[k].
      private double gmax
      The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]).
      private double gmin
      The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]).
      private boolean goodb
      Copy of the goodb parameter.
      private boolean hasConverged
      true if the default convergence criterion is verified.
      private double lqnorm
      The estimate of the norm of P * rL[k-1].
      private RealLinearOperator m
      Reference to the preconditioner, M.
      (package private) static double MACH_PREC
      The machine precision.
      private RealVector mb
      The value of M * b.
      private double minusEpsZeta
      The value of (-eps[k+1] * zeta[k-1]).
      private double oldb
      The value of beta[k].
      private RealVector r1
      The value of beta[k] * M^(-1) * P' * v[k].
      private RealVector r2
      The value of beta[k+1] * M^(-1) * P' * v[k+1].
      private double rnorm
      The value of the updated, preconditioned residual P * r.
      private double shift
      Copy of the shift parameter.
      private double snprod
      The value of s[1] * ...
      private double tnorm
      An estimate of the square of the norm of A * V[k], based on Paige and Saunders (1975), equation (3.3).
      private RealVector wbar
      The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ...
      private RealVector xL
      A reference to the vector to be updated with the solution.
      private RealVector y
      The value of beta[k+1] * P' * v[k+1].
      private double ynorm2
      The value of zeta[1]^2 + ...
    • Field Detail

      • CBRT_MACH_PREC

        static final double CBRT_MACH_PREC
        The cubic root of MACH_PREC.
      • MACH_PREC

        static final double MACH_PREC
        The machine precision.
      • b

        private final RealVector b
        Reference to the right-hand side vector.
      • check

        private final boolean check
        true if symmetry of matrix and conditioner must be checked.
      • delta

        private final double delta
        The value of the custom tolerance δ for the default stopping criterion.
      • beta

        private double beta
        The value of beta[k+1].
      • beta1

        private double beta1
        The value of beta[1].
      • bstep

        private double bstep
        The value of bstep[k-1].
      • cgnorm

        private double cgnorm
        The estimate of the norm of P * rC[k].
      • dbar

        private double dbar
        The value of dbar[k+1] = -beta[k+1] * c[k-1].
      • gammaZeta

        private double gammaZeta
        The value of gamma[k] * zeta[k]. Was called rhs1 in the initial code.
      • gbar

        private double gbar
        The value of gbar[k].
      • gmax

        private double gmax
        The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]).
      • gmin

        private double gmin
        The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]).
      • goodb

        private final boolean goodb
        Copy of the goodb parameter.
      • hasConverged

        private boolean hasConverged
        true if the default convergence criterion is verified.
      • lqnorm

        private double lqnorm
        The estimate of the norm of P * rL[k-1].
      • minusEpsZeta

        private double minusEpsZeta
        The value of (-eps[k+1] * zeta[k-1]). Was called rhs2 in the initial code.
      • mb

        private final RealVector mb
        The value of M * b.
      • oldb

        private double oldb
        The value of beta[k].
      • r1

        private RealVector r1
        The value of beta[k] * M^(-1) * P' * v[k].
      • r2

        private RealVector r2
        The value of beta[k+1] * M^(-1) * P' * v[k+1].
      • rnorm

        private double rnorm
        The value of the updated, preconditioned residual P * r. This value is given by min(cgnorm, lqnorm).
      • shift

        private final double shift
        Copy of the shift parameter.
      • snprod

        private double snprod
        The value of s[1] * ... * s[k-1].
      • tnorm

        private double tnorm
        An estimate of the square of the norm of A * V[k], based on Paige and Saunders (1975), equation (3.3).
      • wbar

        private RealVector wbar
        The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] * v[1]) if goodb is true. Was called w in the initial code.
      • xL

        private final RealVector xL
        A reference to the vector to be updated with the solution. Contains the value of xL[k-1] if goodb is false, (xL[k-1] - bstep[k-1] * v[1]) otherwise.
      • y

        private RealVector y
        The value of beta[k+1] * P' * v[k+1].
      • ynorm2

        private double ynorm2
        The value of zeta[1]^2 + ... + zeta[k-1]^2.
      • bIsNull

        private boolean bIsNull
        The value of b == 0 (exact floating-point equality).
    • Constructor Detail

      • State

        State​(RealLinearOperator a,
              RealLinearOperator m,
              RealVector b,
              boolean goodb,
              double shift,
              double delta,
              boolean check)
        Creates and inits to k = 1 a new instance of this class.
        Parameters:
        a - the linear operator A of the system
        m - the preconditioner, M (can be null)
        b - the right-hand side vector
        goodb - usually false, except if x is expected to contain a large multiple of b
        shift - the amount to be subtracted to all diagonal elements of A
        delta - the δ parameter for the default stopping criterion
        check - true if self-adjointedness of both matrix and preconditioner should be checked
    • Method Detail

      • checkSymmetry

        private static void checkSymmetry​(RealLinearOperator l,
                                          RealVector x,
                                          RealVector y,
                                          RealVector z)
                                   throws NonSelfAdjointOperatorException
        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.
        Parameters:
        l - the linear operator L
        x - the candidate vector x
        y - the candidate vector y = L · x
        z - the vector z = L · y
        Throws:
        NonSelfAdjointOperatorException - when the test fails
      • daxpy

        private static void daxpy​(double a,
                                  RealVector x,
                                  RealVector y)
        A clone of the BLAS DAXPY function, which carries out the operation y ← a · x + y. This is for internal use only: no dimension checks are provided.
        Parameters:
        a - the scalar by which x is to be multiplied
        x - the vector to be added to y
        y - the vector to be incremented
      • daxpbypz

        private static void daxpbypz​(double a,
                                     RealVector x,
                                     double b,
                                     RealVector y,
                                     RealVector z)
        A BLAS-like function, for the operation z ← a · x + b · y + z. This is for internal use only: no dimension checks are provided.
        Parameters:
        a - the scalar by which x is to be multiplied
        x - the first vector to be added to z
        b - the scalar by which y is to be multiplied
        y - the second vector to be added to z
        z - the vector to be incremented
      • refineSolution

        void refineSolution​(RealVector x)

        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.

        Additional upudates are also carried out in case goodb is set to true.

        Parameters:
        x - the vector to be updated with the refined value of xL
      • init

        void init()
        Performs the initial phase of the SYMMLQ algorithm. On return, the value of the state variables of this object correspond to k = 1.
      • update

        void update()
        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 this object correspond to the current iteration count k.
      • updateNorms

        private void updateNorms()
        Computes the norms of the residuals, and checks for convergence. Updates lqnorm and cgnorm.
      • hasConverged

        boolean hasConverged()
        Returns true if the default stopping criterion is fulfilled.
        Returns:
        true if convergence of the iterations has occurred
      • bEqualsNullVector

        boolean bEqualsNullVector()
        Returns true if the right-hand side vector is zero exactly.
        Returns:
        the boolean value of b == 0
      • betaEqualsZero

        boolean betaEqualsZero()
        Returns true if beta is essentially zero. This method is used to check for early stop of the iterations.
        Returns:
        true if beta < MACH_PREC
      • getNormOfResidual

        double getNormOfResidual()
        Returns the norm of the updated, preconditioned residual.
        Returns:
        the norm of the residual, ||P * r||