Forked from
agnumpde / dune-tectonic
1477 commits behind the upstream repository.
-
Elias Pipping authoredElias Pipping authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
test-gradient-method.cc 4.48 KiB
/* -*- mode:c++; mode:semantic -*- */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/exceptions.hh>
#include "samplefunctional.hh"
#include <cassert>
template <int dim, class Function>
double functionTester(
Dune::SampleFunctional<dim, Function> J,
typename Dune::SampleFunctional<dim, Function>::SmallVector &start,
size_t runs) {
typename Dune::SampleFunctional<dim, Function>::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;
}
double const final = J(start);
std::cout << "Final value J(...) = " << final << std::endl;
return final;
}
void testSampleFunction() {
int const dim = 2;
typedef Dune::SampleFunctional<dim, Dune::SampleFunction> SampleFunctional;
SampleFunctional::SmallMatrix A;
A[0][0] = 3;
A[0][1] = 1.5;
A[1][0] = 1.5;
A[1][1] = 4;
SampleFunctional::SmallVector b;
b[0] = 1;
b[1] = 2;
SampleFunctional J(A, b);
SampleFunctional::SmallVector start = b;
start *= 17;
/*
j(x)
= Ax - b + 2/|x| x
= 17*(6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
= (102 - 1 + 2/sqrt(5), 161.5 - 2 + 4/sqrt(5))
*/
SampleFunctional::SmallVector error;
error[0] = -(102 - 1 + 2 / sqrt(5));
error[1] = -(161.5 - 2 + 4 / sqrt(5));
SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error -= returned;
assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
double const ret1 = functionTester(J, start, 6);
start[0] = 279;
start[1] = -96;
double const ret2 = functionTester(J, start, 10);
assert(std::abs(ret1 - ret2) < 1e-5);
}
void testSampleFunctionNonsmooth() {
int const dim = 2;
typedef Dune::SampleFunctional<dim, Dune::SampleFunction> SampleFunctional;
SampleFunctional::SmallMatrix A;
A[0][0] = 3;
A[0][1] = 1.5;
A[1][0] = 1.5;
A[1][1] = 4;
SampleFunctional::SmallVector b;
b[0] = 1;
b[1] = 2;
SampleFunctional J(A, b);
SampleFunctional::SmallVector start;
SampleFunctional::SmallVector error;
/*
for x = b/|b|:
j_(x)
= Ax - b + 1/|x| x
= 1/sqrt(5) (6, 9.5) - (1, 2) + 1/sqrt(5) (1, 2)
= (7/sqrt(5) - 1, 11.5/sqrt(5) - 2)
j+(x)
= Ax - b + 2/|x| x
= 1/sqrt(5) (6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
= (8/sqrt(5) - 1, 13.5/sqrt(5) - 2)
*/
{
start = b;
start /= (b.two_norm() + 1e-12);
assert(start.two_norm() < 1);
SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error[0] = -(7 / sqrt(5) - 1);
error[1] = -(11.5 / sqrt(5) - 2);
error -= returned;
assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
functionTester(J, start, 6);
}
std::cout << std::endl;
{
start = b;
start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
assert(start.two_norm() > 1);
SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error[0] = -(8 / sqrt(5) - 1);
error[1] = -(13.5 / sqrt(5) - 2);
error -= returned;
assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
functionTester(J, start, 6);
}
}
void testTrivialFunction() {
int const dim = 2;
typedef Dune::SampleFunctional<dim> SampleFunctional;
SampleFunctional::SmallMatrix A;
A[0][0] = 3;
A[0][1] = 1.5;
A[1][0] = 1.5;
A[1][1] = 4;
SampleFunctional::SmallVector b;
b[0] = 1;
b[1] = 2;
SampleFunctional J(A, b);
SampleFunctional::SmallVector start = b;
start *= 17;
/*
j(x)
= Ax - b
= 17*(6, 9.5) - (1, 2)
= (102 - 1, 161.5 - 2)
*/
SampleFunctional::SmallVector error;
error[0] = -101;
error[1] = -159.5;
SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error -= returned;
assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
double const ret1 = functionTester(J, start, 6);
std::cout << std::endl;
start[0] = 279;
start[1] = -96;
double const ret2 = functionTester(J, start, 16);
assert(std::abs(ret1 - ret2) < 1e-5);
}
int main() {
try {
testSampleFunction();
std::cout << std::endl << std::endl << std::endl;
testSampleFunctionNonsmooth();
std::cout << std::endl << std::endl << std::endl;
testTrivialFunction();
return 0;
}
catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
}
}