|
|||||||||
PREV CLASS NEXT CLASS | FRAMES NO FRAMES | ||||||||
SUMMARY: NESTED | FIELD | CONSTR | METHOD | DETAIL: FIELD | CONSTR | METHOD |
java.lang.Objectorg.apache.commons.math3.linear.IterativeLinearSolver
org.apache.commons.math3.linear.PreconditionedIterativeLinearSolver
org.apache.commons.math3.linear.SymmLQ
public class SymmLQ
Implementation of the SYMMLQ iterative linear solver proposed by Paige and Saunders (1975). This implementation is largely based on the FORTRAN code by Pr. Michael A. Saunders, available here.
SYMMLQ is designed to solve the system of linear equations A · x = b
where A is an n × n self-adjoint linear operator (defined as a
RealLinearOperator
), and b is a given vector. The operator A is not
required to be positive definite. If A is known to be definite, the method of
conjugate gradients might be preferred, since it will require about the same
number of iterations as SYMMLQ but slightly less work per iteration.
SYMMLQ is designed to solve the system (A - shift · I) · x = b, where shift is a specified scalar value. If shift and b are suitably chosen, the computed vector x may approximate an (unnormalized) eigenvector of A, as in the methods of inverse iteration and/or Rayleigh-quotient iteration. Again, the linear operator (A - shift · I) need not be positive definite (but must be self-adjoint). The work per iteration is very slightly less if shift = 0.
Preconditioning may reduce the number of iterations required. The solver may be provided with a positive definite preconditioner M = C · CT that is known to approximate (A - shift · I) in some sense, where systems of the form M · y = x can be solved efficiently. Then SYMMLQ will implicitly solve the system of equations P · (A - shift · I) · PT · xhat = P · b, i.e. Ahat · xhat = bhat, where P = C-1, Ahat = P · (A - shift · I) · PT, bhat = P · b, and return the solution x = PT · xhat. The associated residual is rhat = bhat - Ahat · xhat = P · [b - (A - shift · I) · x] = P · r.
A default stopping criterion is implemented. The iterations stop when || rhat || ≤ δ || Ahat || || xhat ||, where xhat is the current estimate of the solution of the transformed system, rhat the current estimate of the corresponding residual, and δ a user-specified tolerance.
In the present context, an iteration should be understood as one evaluation of the matrix-vector product A · x. The initialization phase therefore counts as one iteration. If the user requires checks on the symmetry of A, this entails one further matrix-vector product in the initial phase. This further product is not accounted for in the iteration count. In other words, the number of iterations required to reach convergence will be identical, whether checks have been required or not.
The present definition of the iteration count differs from that adopted in the original FOTRAN code, where the initialization phase was not taken into account.
The x
parameter in
solve(RealLinearOperator, RealVector, RealVector)
,solve(RealLinearOperator, RealLinearOperator, RealVector, RealVector)
},solveInPlace(RealLinearOperator, RealVector, RealVector)
,solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector)
,solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector, boolean, double)
,
Besides standard DimensionMismatchException
, this class might throw
NonSelfAdjointOperatorException
if the linear operator or the
preconditioner are not symmetric. In this case, the ExceptionContext
provides more information
"operator"
points to the offending linear operator, say L,"vector1"
points to the first offending vector, say x,
"vector2"
points to the second offending vector, say y, such
that xT · L · y ≠ yT · L
· x (within a certain accuracy).
NonPositiveDefiniteOperatorException
might also be thrown in case the
preconditioner is not positive definite. The relevant keys to the
ExceptionContext
are
"operator"
, which points to the offending linear operator,
say L,"vector"
, which points to the offending vector, say x, such
that xT · L · x < 0.
Nested Class Summary | |
---|---|
private class |
SymmLQ.State
A simple container holding the non-final variables used in the iterations. |
private static class |
SymmLQ.SymmLQEvent
The type of all events fired by this implementation of the SYMMLQ method. |
Field Summary | |
---|---|
private static double |
CBRT_MACH_PREC
The cubic root of MACH_PREC . |
private boolean |
check
true if symmetry of matrix and conditioner must be checked. |
private double |
delta
The value of the custom tolerance δ for the default stopping criterion. |
private static double |
MACH_PREC
The machine precision. |
private static String |
OPERATOR
Key for the exception context. |
private static String |
THRESHOLD
Key for the exception context. |
private static String |
VECTOR
Key for the exception context. |
private static String |
VECTOR1
Key for the exception context. |
private static String |
VECTOR2
Key for the exception context. |
Constructor Summary | |
---|---|
SymmLQ(int maxIterations,
double delta,
boolean check)
Creates a new instance of this class, with default stopping criterion. |
|
SymmLQ(IterationManager manager,
double delta,
boolean check)
Creates a new instance of this class, with default stopping criterion and custom iteration manager. |
Method Summary | |
---|---|
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 BLAS DAXPY function, which carries out the
operation y ← a · x + y. |
boolean |
getCheck()
Returns true if symmetry of the matrix, and symmetry as well as
positive definiteness of the preconditioner should be checked. |
RealVector |
solve(RealLinearOperator a,
RealLinearOperator minv,
RealVector b)
Returns an estimate of the solution to the linear system A · x = b. |
RealVector |
solve(RealLinearOperator a,
RealLinearOperator minv,
RealVector b,
boolean goodb,
double shift)
Returns an estimate of the solution to the linear system (A - shift · I) · x = b. |
RealVector |
solve(RealLinearOperator a,
RealLinearOperator minv,
RealVector b,
RealVector x)
Returns an estimate of the solution to the linear system A · x = b. |
RealVector |
solve(RealLinearOperator a,
RealVector b)
Returns an estimate of the solution to the linear system A · x = b. |
RealVector |
solve(RealLinearOperator a,
RealVector b,
boolean goodb,
double shift)
Returns the solution to the system (A - shift · I) · x = b. |
RealVector |
solve(RealLinearOperator a,
RealVector b,
RealVector x)
Returns an estimate of the solution to the linear system A · x = b. |
RealVector |
solveInPlace(RealLinearOperator a,
RealLinearOperator minv,
RealVector b,
RealVector x)
Returns an estimate of the solution to the linear system A · x = b. |
RealVector |
solveInPlace(RealLinearOperator a,
RealLinearOperator minv,
RealVector b,
RealVector x,
boolean goodb,
double shift)
Returns an estimate of the solution to the linear system (A - shift · I) · x = b. |
RealVector |
solveInPlace(RealLinearOperator a,
RealVector b,
RealVector x)
Returns an estimate of the solution to the linear system A · x = b. |
private static void |
throwNPDLOException(RealLinearOperator l,
RealVector v)
Throws a new NonPositiveDefiniteOperatorException with
appropriate context. |
Methods inherited from class org.apache.commons.math3.linear.PreconditionedIterativeLinearSolver |
---|
checkParameters |
Methods inherited from class org.apache.commons.math3.linear.IterativeLinearSolver |
---|
checkParameters, getIterationManager |
Methods inherited from class java.lang.Object |
---|
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait |
Field Detail |
---|
private static final double CBRT_MACH_PREC
MACH_PREC
.
private static final double MACH_PREC
private static final String OPERATOR
private static final String THRESHOLD
private static final String VECTOR
private static final String VECTOR1
private static final String VECTOR2
private final boolean check
true
if symmetry of matrix and conditioner must be checked.
private final double delta
Constructor Detail |
---|
public SymmLQ(int maxIterations, double delta, boolean check)
check
to true
entails an extra matrix-vector product in the initial phase.
maxIterations
- the maximum number of iterationsdelta
- the δ parameter for the default stopping criterioncheck
- true
if self-adjointedness of both matrix and
preconditioner should be checkedpublic SymmLQ(IterationManager manager, double delta, boolean check)
check
to true
entails an extra matrix-vector product in
the initial phase.
manager
- the custom iteration managerdelta
- the δ parameter for the default stopping criterioncheck
- true
if self-adjointedness of both matrix and
preconditioner should be checkedMethod Detail |
---|
private static void checkSymmetry(RealLinearOperator l, RealVector x, RealVector y, RealVector z) throws NonSelfAdjointOperatorException
l
- the linear operator Lx
- the candidate vector xy
- the candidate vector y = L · xz
- the vector z = L · y
NonSelfAdjointOperatorException
- when the test failsprivate static void daxpbypz(double a, RealVector x, double b, RealVector y, RealVector z)
a
- the scalar by which x
is to be multipliedx
- the first vector to be added to z
b
- the scalar by which y
is to be multipliedy
- the second vector to be added to z
z
- the vector to be incrementedprivate static void daxpy(double a, RealVector x, RealVector y)
DAXPY
function, which carries out the
operation y ← a · x + y. This is for internal use only: no
dimension checks are provided.
a
- the scalar by which x
is to be multipliedx
- the vector to be added to y
y
- the vector to be incrementedprivate static void throwNPDLOException(RealLinearOperator l, RealVector v) throws NonPositiveDefiniteOperatorException
NonPositiveDefiniteOperatorException
with
appropriate context.
l
- the offending linear operatorv
- the offending vector
NonPositiveDefiniteOperatorException
- in any circumstancespublic final boolean getCheck()
true
if symmetry of the matrix, and symmetry as well as
positive definiteness of the preconditioner should be checked.
true
if the tests are to be performedpublic RealVector solve(RealLinearOperator a, RealLinearOperator minv, RealVector b) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, MaxCountExceededException, NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException, IllConditionedOperatorException
solve
in class PreconditionedIterativeLinearSolver
a
- the linear operator A of the systemminv
- the inverse of the preconditioner, M-1
(can be null
)b
- the right-hand side vector
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
or minv
is not self-adjoint
NonPositiveDefiniteOperatorException
- if minv
is not
positive definite
IllConditionedOperatorException
- if a
is ill-conditioned
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
or minv
is not
square
DimensionMismatchException
- if minv
or b
have
dimensions inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at constructionpublic RealVector solve(RealLinearOperator a, RealLinearOperator minv, RealVector b, boolean goodb, double shift) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, MaxCountExceededException, NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException, IllConditionedOperatorException
If the solution x is expected to contain a large multiple of b
(as in Rayleigh-quotient iteration), then better precision may be
achieved with goodb
set to true
; this however requires an
extra call to the preconditioner.
shift
should be zero if the system A · x = b is to be
solved. Otherwise, it could be an approximation to an eigenvalue of A,
such as the Rayleigh quotient bT · A · b /
(bT · b) corresponding to the vector b. If b is
sufficiently like an eigenvector corresponding to an eigenvalue near
shift, then the computed x may have very large components. When
normalized, x may be closer to an eigenvector than b.
a
- the linear operator A of the systemminv
- the inverse of the preconditioner, M-1
(can be null
)b
- the right-hand side vectorgoodb
- 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
x
(shallow copy)
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
or minv
is not
square
DimensionMismatchException
- if minv
or b
have
dimensions inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at construction
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
or minv
is not self-adjoint
NonPositiveDefiniteOperatorException
- if minv
is not
positive definite
IllConditionedOperatorException
- if a
is ill-conditionedpublic RealVector solve(RealLinearOperator a, RealLinearOperator minv, RealVector b, RealVector x) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException, IllConditionedOperatorException, MaxCountExceededException
solve
in class PreconditionedIterativeLinearSolver
x
- not meaningful in this implementation; should not be considered
as an initial guess (more)a
- the linear operator A of the systemminv
- the inverse of the preconditioner, M-1
(can be null
)b
- the right-hand side vector
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
or minv
is not self-adjoint
NonPositiveDefiniteOperatorException
- if minv
is not
positive definite
IllConditionedOperatorException
- if a
is ill-conditioned
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
or minv
is not
square
DimensionMismatchException
- if minv
, b
or
x0
have dimensions inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at constructionpublic RealVector solve(RealLinearOperator a, RealVector b) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, IllConditionedOperatorException, MaxCountExceededException
solve
in class PreconditionedIterativeLinearSolver
a
- the linear operator A of the systemb
- the right-hand side vector
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
is not self-adjoint
IllConditionedOperatorException
- if a
is ill-conditioned
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
is not square
DimensionMismatchException
- if b
has dimensions
inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at constructionpublic RealVector solve(RealLinearOperator a, RealVector b, boolean goodb, double shift) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, IllConditionedOperatorException, MaxCountExceededException
If the solution x is expected to contain a large multiple of b
(as in Rayleigh-quotient iteration), then better precision may be
achieved with goodb
set to true
.
shift
should be zero if the system A · x = b is to be
solved. Otherwise, it could be an approximation to an eigenvalue of A,
such as the Rayleigh quotient bT · A · b /
(bT · b) corresponding to the vector b. If b is
sufficiently like an eigenvector corresponding to an eigenvalue near
shift, then the computed x may have very large components. When
normalized, x may be closer to an eigenvector than b.
a
- the linear operator A of the systemb
- the right-hand side vectorgoodb
- 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
x
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
is not square
DimensionMismatchException
- if b
has dimensions
inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at construction
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
is not self-adjoint
IllConditionedOperatorException
- if a
is ill-conditionedpublic RealVector solve(RealLinearOperator a, RealVector b, RealVector x) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, IllConditionedOperatorException, MaxCountExceededException
solve
in class PreconditionedIterativeLinearSolver
x
- not meaningful in this implementation; should not be considered
as an initial guess (more)a
- the linear operator A of the systemb
- the right-hand side vector
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
is not self-adjoint
IllConditionedOperatorException
- if a
is ill-conditioned
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
is not square
DimensionMismatchException
- if b
or x0
have
dimensions inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at constructionpublic RealVector solveInPlace(RealLinearOperator a, RealLinearOperator minv, RealVector b, RealVector x) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException, IllConditionedOperatorException, MaxCountExceededException
solveInPlace
in class PreconditionedIterativeLinearSolver
x
- the vector to be updated with the solution; x
should
not be considered as an initial guess (more)a
- the linear operator A of the systemminv
- the inverse of the preconditioner, M-1
(can be null
)b
- the right-hand side vector
x0
(shallow copy) updated with the
solution
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
or minv
is not self-adjoint
NonPositiveDefiniteOperatorException
- if minv
is not
positive definite
IllConditionedOperatorException
- if a
is ill-conditioned
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
or minv
is not
square
DimensionMismatchException
- if minv
, b
or
x0
have dimensions inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at construction.public RealVector solveInPlace(RealLinearOperator a, RealLinearOperator minv, RealVector b, RealVector x, boolean goodb, double shift) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, NonPositiveDefiniteOperatorException, IllConditionedOperatorException, MaxCountExceededException
If the solution x is expected to contain a large multiple of b
(as in Rayleigh-quotient iteration), then better precision may be
achieved with goodb
set to true
; this however requires an
extra call to the preconditioner.
shift
should be zero if the system A · x = b is to be
solved. Otherwise, it could be an approximation to an eigenvalue of A,
such as the Rayleigh quotient bT · A · b /
(bT · b) corresponding to the vector b. If b is
sufficiently like an eigenvector corresponding to an eigenvalue near
shift, then the computed x may have very large components. When
normalized, x may be closer to an eigenvector than b.
a
- the linear operator A of the systemminv
- the inverse of the preconditioner, M-1
(can be null
)b
- the right-hand side vectorx
- the vector to be updated with the solution; x
should
not be considered as an initial guess (more)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
x
(shallow copy).
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
or minv
is not
square
DimensionMismatchException
- if minv
, b
or
x
have dimensions inconsistent with a
.
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at construction
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
or minv
is not self-adjoint
NonPositiveDefiniteOperatorException
- if minv
is not
positive definite
IllConditionedOperatorException
- if a
is ill-conditionedpublic RealVector solveInPlace(RealLinearOperator a, RealVector b, RealVector x) throws NullArgumentException, NonSquareOperatorException, DimensionMismatchException, NonSelfAdjointOperatorException, IllConditionedOperatorException, MaxCountExceededException
solveInPlace
in class PreconditionedIterativeLinearSolver
x
- the vector to be updated with the solution; x
should
not be considered as an initial guess (more)a
- the linear operator A of the systemb
- the right-hand side vector
x0
(shallow copy) updated with the
solution
NonSelfAdjointOperatorException
- if getCheck()
is
true
, and a
is not self-adjoint
IllConditionedOperatorException
- if a
is ill-conditioned
NullArgumentException
- if one of the parameters is null
NonSquareOperatorException
- if a
is not square
DimensionMismatchException
- if b
or x0
have
dimensions inconsistent with a
MaxCountExceededException
- at exhaustion of the iteration count,
unless a custom
callback
has been set at construction
|
|||||||||
PREV CLASS NEXT CLASS | FRAMES NO FRAMES | ||||||||
SUMMARY: NESTED | FIELD | CONSTR | METHOD | DETAIL: FIELD | CONSTR | METHOD |