ConjugateGradient extends PreconditionedIterativeLinearSolver { @Override public RealVector solveInPlace(final RealLinearOperator a, final RealLinearOperator m, final RealVector b, final RealVector x0) throws NullArgumentException, NonPositiveDefiniteOperatorException, NonSquareOperatorException, DimensionMismatchException, MaxCountExceededException { checkParameters(a, m, b, x0); final IterationManager manager = getIterationManager(); manager.resetIterationCount(); final double rmax = delta * b.getNorm(); final RealVector bro = RealVector.unmodifiableRealVector(b); manager.incrementIterationCount(); 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 (m == 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 (m != null) { z = m.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, m); 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; } } } ConjugateGradient(final int maxIterations, final double delta,
                             final boolean check);  ConjugateGradient(final IterationManager manager,
                             final double delta, final boolean check); final boolean getCheck(); @Override RealVector solveInPlace(final RealLinearOperator a,
                                   final RealLinearOperator m,
                                   final RealVector b,
                                   final RealVector x0); static final String OPERATOR; static final String VECTOR; }
@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 ConjugateGradient(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-10 * FastMath.abs(expected); final String msg = String.format("entry[%d][%d)", i, j); Assert.assertEquals(msg, expected, actual, delta); } } }