public RealVector solveInPlace(final RealLinearOperator a,
final RealLinearOperator minv, final RealVector b, final RealVector x0)
throws NullArgumentException, NonSquareOperatorException,
DimensionMismatchException, MaxCountExceededException {
checkParameters(a, minv, b, x0);
final IterationManager manager = getIterationManager();
// Initialization of default stopping criterion
manager.resetIterationCount();
final double rmax = delta * b.getNorm();
final RealVector bro = RealVector.unmodifiableRealVector(b);
// Initialization phase counts as one iteration.
manager.incrementIterationCount();
// p and x are constructed as copies of x0, since presumably, the type
// of x is optimized for the calculation of the matrix-vector product
// A.x.
final RealVector x = x0;
final RealVector xro = RealVector.unmodifiableRealVector(x);
final RealVector p = x.copy();
RealVector q = a.operate(p);
final RealVector r = b.combine(1, -1, q);
final RealVector rro = RealVector.unmodifiableRealVector(r);
double rnorm = r.getNorm();
RealVector z;
if (minv == null) {
z = r;
} else {
z = null;
}
IterativeLinearSolverEvent evt;
evt = new DefaultIterativeLinearSolverEvent(this,
manager.getIterations(), xro, bro, rro, rnorm);
manager.fireInitializationEvent(evt);
if (rnorm <= rmax) {
manager.fireTerminationEvent(evt);
return x;
}
double rhoPrev = 0.;
while (true) {
manager.incrementIterationCount();
evt = new DefaultIterativeLinearSolverEvent(this,
manager.getIterations(), xro, bro, rro, rnorm);
manager.fireIterationStartedEvent(evt);
if (minv != null) {
z = minv.operate(r);
}
final double rhoNext = r.dotProduct(z);
if (check && (rhoNext <= 0.)) {
final NonPositiveDefiniteOperatorException e;
e = new NonPositiveDefiniteOperatorException();
final ExceptionContext context = e.getContext();
context.setValue(OPERATOR, minv);
context.setValue(VECTOR, r);
throw e;
}
if (manager.getIterations() == 2) {
p.setSubVector(0, z);
} else {
p.combineToSelf(rhoNext / rhoPrev, 1., z);
}
q = a.operate(p);
final double pq = p.dotProduct(q);
if (check && (pq <= 0.)) {
final NonPositiveDefiniteOperatorException e;
e = new NonPositiveDefiniteOperatorException();
final ExceptionContext context = e.getContext();
context.setValue(OPERATOR, a);
context.setValue(VECTOR, p);
throw e;
}
final double alpha = rhoNext / pq;
x.combineToSelf(1., alpha, p);
r.combineToSelf(1., -alpha, q);
rhoPrev = rhoNext;
rnorm = r.getNorm();
evt = new DefaultIterativeLinearSolverEvent(this,
manager.getIterations(), xro, bro, rro, rnorm);
manager.fireIterationPerformedEvent(evt);
if (rnorm <= rmax) {
manager.fireTerminationEvent(evt);
return x;
}
}
}