diff --git a/src/myblockproblem.hh b/src/myblockproblem.hh index acf4801fba56ea33b7e7d51e8315c1ac82164e0d..b8c2415d235072c2c8e76e785b56501e5153e37c 100644 --- a/src/myblockproblem.hh +++ b/src/myblockproblem.hh @@ -32,6 +32,7 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem { class IterateObject; MyBlockProblem(MyConvexProblemType& problem) : problem(problem) { + // TODO: Is it clever to create a bisection here? bisection = Bisection(0.0, 1.0, 1e-15, true, 1e-14); }; @@ -118,10 +119,7 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject { Dune::SampleFunctional<block_size> localJ(*localA, localb, phi); LocalVectorType correction; - for (size_t i = 1; i <= 10; ++i) { // FIXME: hardcoded value - Dune::minimise(localJ, ui, correction); - ui += correction; - } + Dune::minimise(localJ, ui, 10); // FIXME: hardcoded value } } diff --git a/src/samplefunctional.hh b/src/samplefunctional.hh index 0a1dc654abee63e58c9b361e3306b00801cd91d1..9b7a710047b7b932635da07e6956b89abfba1cb0 100644 --- a/src/samplefunctional.hh +++ b/src/samplefunctional.hh @@ -118,8 +118,8 @@ template <int dimension> class SampleFunctional { }; template <class Functional> -void minimise(const Functional J, const typename Functional::SmallVector x, - typename Functional::SmallVector &corr, +void minimise(const Functional J, typename Functional::SmallVector &x, + size_t steps = 1, Bisection const &bisection = Bisection(0.0, // acceptError: Stop if the search interval has // become smaller than this number @@ -130,59 +130,60 @@ void minimise(const Functional J, const typename Functional::SmallVector x, // minimization { typedef typename Functional::SmallVector SmallVector; - SmallVector descDir; - J.descentDirection(x, descDir); - if (descDir == SmallVector(0.0)) { - corr = SmallVector(0.0); - return; - } + for (size_t step = 0; step < steps; ++step) { + SmallVector descDir; + J.descentDirection(x, descDir); - // {{{ Construct a restriction of J to the line x + t * descDir + if (descDir == SmallVector(0.0)) + return; - /* We have + // {{{ Construct a restriction of J to the line x + t * descDir - 1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u> + /* We have - since A is symmetric. - */ - SmallVector tmp; + 1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u> - J.A.mv(descDir, tmp); // Av - double const JRestA = tmp * descDir; // <Av,v> + since A is symmetric. + */ + SmallVector tmp; - J.A.mv(x, tmp); // Au - double const JRestb = (J.b - tmp) * descDir; // <b-Au,v> + J.A.mv(descDir, tmp); // Av + double const JRestA = tmp * descDir; // <Av,v> - typedef typename Functional::NonlinearityType MyNonlinearityType; - MyNonlinearityType phi = J.phi; - typedef DirectionalConvexFunction<MyNonlinearityType> - MyDirectionalConvexFunctionType; - // FIXME: We cannot pass J.phi directly because the constructor - // does not allow for constant arguments - MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir); - // }}} + J.A.mv(x, tmp); // Au + double const JRestb = (J.b - tmp) * descDir; // <b-Au,v> - { // Debug - Interval<double> D; - JRest.subDiff(0, D); + typedef typename Functional::NonlinearityType MyNonlinearityType; + MyNonlinearityType phi = J.phi; + typedef DirectionalConvexFunction<MyNonlinearityType> + MyDirectionalConvexFunctionType; + // FIXME: We cannot pass J.phi directly because the constructor + // does not allow for constant arguments + MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir); + // }}} - dverb - << "## Directional derivative (as per subdifferential of restriction): " - << D[1] << " (coordinates of the restriction)" << std::endl; - assert(D[1] <= - 0); // We should not be minimising in this direction otherwise - } + { // Debug + Interval<double> D; + JRest.subDiff(0, D); - int count; - // FIXME: The value of x_old should not matter if the factor is 1.0, correct? - double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count); - dverb << "Number of iterations in the bisection method: " << count - << std::endl; - ; + dverb << "## Directional derivative (as per subdifferential of " + "restriction): " << D[1] << " (coordinates of the restriction)" + << std::endl; + assert(D[1] <= + 0); // We should not be minimising in this direction otherwise + } + + int count; + // FIXME: The value of x_old should not matter if the factor is 1.0, + // correct? + double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count); + dverb << "Number of iterations in the bisection method: " << count + << std::endl; + ; - corr = descDir; - corr *= stepsize; + x.axpy(stepsize, descDir); + } } } #endif diff --git a/src/test-gradient-method.cc b/src/test-gradient-method.cc index 8165981f961ef9b0e73451c2d3553e126c472850..8a1b6dcf83e63d50f8593c6ebb90fb6f2417c265 100644 --- a/src/test-gradient-method.cc +++ b/src/test-gradient-method.cc @@ -15,14 +15,8 @@ template <int dim> double functionTester(Dune::SampleFunctional<dim> J, typename Dune::SampleFunctional<dim>::SmallVector &start, size_t runs) { - typename Dune::SampleFunctional<dim>::SmallVector correction; std::cout << "Old value: J(...) = " << J(start) << std::endl; - for (size_t i = 1; i <= runs; ++i) { - Dune::minimise(J, start, correction); - start += correction; - if (i != runs) - std::cout << "New value: J(...) = " << J(start) << std::endl; - } + Dune::minimise(J, start, runs); double const final = J(start); std::cout << "Final value J(...) = " << final << std::endl; return final;