Skip to content
Snippets Groups Projects
Commit b2c65d4f authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Kill tangentialMinimisation

parent 9cfd8843
No related branches found
No related tags found
No related merge requests found
#ifndef CIRCULAR_CONVEX_FUNCTION_HH
#define CIRCULAR_CONVEX_FUNCTION_HH
#include <cmath>
#include <dune/fufem/interval.hh>
namespace Dune {
template <class MatrixType, class VectorType> class CircularConvexFunction {
public:
CircularConvexFunction(MatrixType const &A, VectorType const &b,
VectorType const &x, VectorType const &dir)
: A(A),
b(b),
x(x),
dir(dir),
xnorm(x.two_norm()),
dnorm(dir.two_norm()) {}
double quadraticPart() const { return 0; }
double linearPart() const { return 0; }
void subDiff(double m, Interval<double> &D) const {
VectorType x;
cartesian(m, x);
VectorType t;
tangent(m, t);
VectorType tmp;
A.mv(x, tmp); // Ax
tmp -= b; // Ax - b
D[0] = D[1] = tmp * t; // <Ax - b,t>
}
void domain(Interval<double> &domain) const {
domain[0] = -2 * M_PI;
domain[1] = 2 * M_PI;
}
void cartesian(double m, VectorType &y) const {
y = 0;
y.axpy(std::cos(m), x);
y.axpy(std::sin(m) * xnorm / dnorm, dir);
}
private:
MatrixType const &A;
VectorType const &b;
VectorType const &x;
VectorType const &dir;
double const dnorm;
double const xnorm;
/* If x and d were normalised, we would have
cartesian(a) = cos(a) * x + sin(a) * d and
tangent(a) = -sin(a) * x + cos(a) * d.
Since we x and d are not normalised and the return of
cartesian() is fixed, we scale the tangent.
*/
void tangent(double m, VectorType &y) const {
y = 0;
y.axpy(-std::sin(m) * dnorm, x);
y.axpy(std::cos(m) * xnorm, dir);
}
};
}
#endif
...@@ -8,7 +8,6 @@ ...@@ -8,7 +8,6 @@
#include <dune/fufem/interval.hh> #include <dune/fufem/interval.hh>
#include <dune/tnnmg/problem-classes/bisection.hh> #include <dune/tnnmg/problem-classes/bisection.hh>
#include "circularconvexfunction.hh"
#include "mydirectionalconvexfunction.hh" #include "mydirectionalconvexfunction.hh"
namespace Dune { namespace Dune {
...@@ -70,31 +69,6 @@ void descentMinimisation(Functional const &J, ...@@ -70,31 +69,6 @@ void descentMinimisation(Functional const &J,
x.axpy(stepsize, descDir); x.axpy(stepsize, descDir);
} }
template <class Functional>
void tangentialMinimisation(Functional const &J,
typename Functional::SmallVector &x,
typename Functional::SmallVector const &descDir,
Bisection const &bisection) {
using SmallMatrix = typename Functional::SmallMatrix;
using SmallVector = typename Functional::SmallVector;
// We completely ignore the nonlinearity here -- when restricted
// to a circle, it just enters as a constant!
CircularConvexFunction<SmallMatrix, SmallVector> const JRest(J.A, J.b, x,
descDir);
int count;
double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
dverb << "Number of iterations in the bisection method: " << count
<< std::endl;
;
// Since x is used in the computation of the rhs, do not write to it directly
SmallVector tmp;
JRest.cartesian(stepsize, tmp);
x = tmp;
}
template <class Functional> template <class Functional>
void minimise(Functional const &J, typename Functional::SmallVector &x, void minimise(Functional const &J, typename Functional::SmallVector &x,
size_t steps, Bisection const &bisection) { size_t steps, Bisection const &bisection) {
...@@ -102,19 +76,12 @@ void minimise(Functional const &J, typename Functional::SmallVector &x, ...@@ -102,19 +76,12 @@ void minimise(Functional const &J, typename Functional::SmallVector &x,
for (size_t step = 0; step < steps; ++step) { for (size_t step = 0; step < steps; ++step) {
SmallVector descDir; SmallVector descDir;
bool const linesearchp = J.descentDirection(x, descDir); J.descentDirection(x, descDir);
if (descDir.two_norm() < 1e-14) // TODO: Make controllable if (descDir.two_norm() < 1e-14) // TODO: Make controllable
return; return;
if (linesearchp) { descentMinimisation(J, x, descDir, bisection);
descentMinimisation(J, x, descDir, bisection);
} else {
Bisection slowBisection(bisection);
slowBisection.setFastQuadratic(false);
tangentialMinimisation(J, x, descDir, slowBisection);
}
} }
} }
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment