Class SymmLQ.State
- Enclosing class:
SymmLQ
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
FieldsModifier and TypeFieldDescriptionprivate final RealLinearOperator
Reference to the linear operator.private final RealVector
Reference to the right-hand side vector.private double
The value of beta[k+1].private double
The value of beta[1].private boolean
The value ofb == 0
(exact floating-point equality).private double
The value of bstep[k-1].(package private) static final double
The cubic root ofMACH_PREC
.private double
The estimate of the norm of P * rC[k].private final boolean
true
if symmetry of matrix and conditioner must be checked.private double
The value of dbar[k+1] = -beta[k+1] * c[k-1].private final double
The value of the custom tolerance δ for the default stopping criterion.private double
The value of gamma[k] * zeta[k].private double
The value of gbar[k].private double
The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]).private double
The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]).private final boolean
Copy of thegoodb
parameter.private boolean
true
if the default convergence criterion is verified.private double
The estimate of the norm of P * rL[k-1].private final RealLinearOperator
Reference to the preconditioner, M.(package private) static final double
The machine precision.private final RealVector
The value of M * b.private double
The value of (-eps[k+1] * zeta[k-1]).private double
The value of beta[k].private RealVector
The value of beta[k] * M^(-1) * P' * v[k].private RealVector
The value of beta[k+1] * M^(-1) * P' * v[k+1].private double
The value of the updated, preconditioned residual P * r.private final double
Copy of theshift
parameter.private double
The value of s[1] * ...private double
An estimate of the square of the norm of A * V[k], based on Paige and Saunders (1975), equation (3.3).private RealVector
The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ...private final RealVector
A reference to the vector to be updated with the solution.private RealVector
The value of beta[k+1] * P' * v[k+1].private double
The value of zeta[1]^2 + ... -
Constructor Summary
ConstructorsConstructorDescriptionState
(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. -
Method Summary
Modifier and TypeMethodDescription(package private) boolean
Returnstrue
if the right-hand side vector is zero exactly.(package private) boolean
Returnstrue
ifbeta
is essentially zero.private static void
checkSymmetry
(RealLinearOperator l, RealVector x, RealVector y, RealVector z) Performs a symmetry check on the specified linear operator, and throws an exception in case this check fails.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.private static void
daxpy
(double a, RealVector x, RealVector y) A clone of the BLASDAXPY
function, which carries out the operation y ← a · x + y.(package private) double
Returns the norm of the updated, preconditioned residual.(package private) boolean
Returnstrue
if the default stopping criterion is fulfilled.(package private) void
init()
Performs the initial phase of the SYMMLQ algorithm.(package private) void
Move to the CG point if it seems better.private static void
Throws a newNonPositiveDefiniteOperatorException
with appropriate context.(package private) void
update()
Performs the next iteration of the algorithm.private void
Computes the norms of the residuals, and checks for convergence.
-
Field Details
-
CBRT_MACH_PREC
static final double CBRT_MACH_PRECThe cubic root ofMACH_PREC
. -
MACH_PREC
static final double MACH_PRECThe machine precision. -
a
Reference to the linear operator. -
b
Reference to the right-hand side vector. -
check
private final boolean checktrue
if symmetry of matrix and conditioner must be checked. -
delta
private final double deltaThe value of the custom tolerance δ for the default stopping criterion. -
beta
private double betaThe value of beta[k+1]. -
beta1
private double beta1The value of beta[1]. -
bstep
private double bstepThe value of bstep[k-1]. -
cgnorm
private double cgnormThe estimate of the norm of P * rC[k]. -
dbar
private double dbarThe value of dbar[k+1] = -beta[k+1] * c[k-1]. -
gammaZeta
private double gammaZetaThe value of gamma[k] * zeta[k]. Was calledrhs1
in the initial code. -
gbar
private double gbarThe value of gbar[k]. -
gmax
private double gmaxThe value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). -
gmin
private double gminThe value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). -
goodb
private final boolean goodbCopy of thegoodb
parameter. -
hasConverged
private boolean hasConvergedtrue
if the default convergence criterion is verified. -
lqnorm
private double lqnormThe estimate of the norm of P * rL[k-1]. -
m
Reference to the preconditioner, M. -
minusEpsZeta
private double minusEpsZetaThe value of (-eps[k+1] * zeta[k-1]). Was calledrhs2
in the initial code. -
mb
The value of M * b. -
oldb
private double oldbThe value of beta[k]. -
r1
The value of beta[k] * M^(-1) * P' * v[k]. -
r2
The value of beta[k+1] * M^(-1) * P' * v[k+1]. -
rnorm
private double rnorm -
shift
private final double shiftCopy of theshift
parameter. -
snprod
private double snprodThe value of s[1] * ... * s[k-1]. -
tnorm
private double tnormAn estimate of the square of the norm of A * V[k], based on Paige and Saunders (1975), equation (3.3). -
wbar
The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] * v[1]) ifgoodb
istrue
. Was calledw
in the initial code. -
xL
A reference to the vector to be updated with the solution. Contains the value of xL[k-1] ifgoodb
isfalse
, (xL[k-1] - bstep[k-1] * v[1]) otherwise. -
y
The value of beta[k+1] * P' * v[k+1]. -
ynorm2
private double ynorm2The value of zeta[1]^2 + ... + zeta[k-1]^2. -
bIsNull
private boolean bIsNullThe value ofb == 0
(exact floating-point equality).
-
-
Constructor Details
-
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 systemm
- the preconditioner, M (can benull
)b
- the right-hand side vectorgoodb
- usuallyfalse
, except ifx
is expected to contain a large multiple ofb
shift
- the amount to be subtracted to all diagonal elements of Adelta
- the δ parameter for the default stopping criterioncheck
-true
if self-adjointedness of both matrix and preconditioner should be checked
-
-
Method Details
-
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 Lx
- the candidate vector xy
- the candidate vector y = L · xz
- the vector z = L · y- Throws:
NonSelfAdjointOperatorException
- when the test fails
-
throwNPDLOException
private static void throwNPDLOException(RealLinearOperator l, RealVector v) throws NonPositiveDefiniteOperatorException Throws a newNonPositiveDefiniteOperatorException
with appropriate context.- Parameters:
l
- the offending linear operatorv
- the offending vector- Throws:
NonPositiveDefiniteOperatorException
- in any circumstances
-
daxpy
A clone of the BLASDAXPY
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 whichx
is to be multipliedx
- the vector to be added toy
y
- the vector to be incremented
-
daxpbypz
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 whichx
is to be multipliedx
- the first vector to be added toz
b
- the scalar by whichy
is to be multipliedy
- the second vector to be added toz
z
- the vector to be incremented
-
refineSolution
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 totrue
.- 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 ofthis
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 ofthis
object correspond to the current iteration countk
. -
updateNorms
private void updateNorms() -
hasConverged
boolean hasConverged()Returnstrue
if the default stopping criterion is fulfilled.- Returns:
true
if convergence of the iterations has occurred
-
bEqualsNullVector
boolean bEqualsNullVector()Returnstrue
if the right-hand side vector is zero exactly.- Returns:
- the boolean value of
b == 0
-
betaEqualsZero
boolean betaEqualsZero()Returnstrue
ifbeta
is essentially zero. This method is used to check for early stop of the iterations.- Returns:
true
ifbeta <
MACH_PREC
-
getNormOfResidual
double getNormOfResidual()Returns the norm of the updated, preconditioned residual.- Returns:
- the norm of the residual, ||P * r||
-