From d457f3e031d889e345b89ddfcbc7327cd10b6895 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 11 Mar 2012 16:31:09 +0100
Subject: [PATCH] Test LambertW function from the GSL
---
src/Makefile.am | 2 ++
src/compute_state.cc | 8 ++++++++
src/compute_state.hh | 2 ++
src/one-body-sample.cc | 20 ++++++++++++++++++--
4 files changed, 30 insertions(+), 2 deletions(-)
diff --git a/src/Makefile.am b/src/Makefile.am
index 85f5ef83..6ec31eef 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -26,6 +26,8 @@ one_body_sample_SOURCES = \
one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc
one_body_sample_CPPFLAGS = \
$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\"
+one_body_sample_LDADD = \
+ $(LDADD) -lgsl
test_gradient_method_SOURCES = \
test-gradient-method.cc
diff --git a/src/compute_state.cc b/src/compute_state.cc
index 15c13a11..46163f4c 100644
--- a/src/compute_state.cc
+++ b/src/compute_state.cc
@@ -1,5 +1,7 @@
#include "LambertW.h"
+#include <gsl/gsl_sf_lambert.h>
+
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/fufem/interval.hh>
@@ -43,3 +45,9 @@ double compute_state_update_lambert(double h, double unorm, double L,
double const rhs = unorm / L - old_state;
return LambertW(0, h * exp(rhs)) - rhs;
}
+
+double compute_state_update_lambert_gsl(double h, double unorm, double L,
+ double old_state) {
+ double const rhs = unorm / L - old_state;
+ return gsl_sf_lambert_W0(h * std::exp(rhs)) - rhs;
+}
diff --git a/src/compute_state.hh b/src/compute_state.hh
index 4ed5680e..2710f625 100644
--- a/src/compute_state.hh
+++ b/src/compute_state.hh
@@ -5,5 +5,7 @@ double compute_state_update_bisection(double h, double unorm, double L,
double old_state);
double compute_state_update_lambert(double h, double unorm, double L,
double old_state);
+double compute_state_update_lambert_gsl(double h, double unorm, double L,
+ double old_state);
#endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 124786b2..6bdedb59 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -360,7 +360,14 @@ int main(int argc, char *argv[]) {
assert(std::abs(1.0 / h * ret2 - (s4_old[i] - unorm / L) / h -
exp(-ret2)) < 1e-10);
- assert(std::abs(ret2 - ret1) < 1e-14);
+ double ret3 =
+ compute_state_update_lambert_gsl(h, unorm, L, s4_old[i][0]);
+ assert(std::abs(1.0 / h * ret3 - (s4_old[i] - unorm / L) / h -
+ exp(-ret3)) < 1e-10);
+
+ assert(std::abs(ret1 - ret2) < 1e-14);
+ assert(std::abs(ret1 - ret3) < 1e-14);
+ assert(std::abs(ret2 - ret3) < 1e-14);
(*s4_new)[i][0] = ret1;
}
@@ -420,11 +427,20 @@ int main(int argc, char *argv[]) {
compute_state_update_bisection(h, unorm, L, s5_old[i][0]);
assert(std::abs(1.0 / h * ret1 - (s5_old[i] - unorm / L) / h -
exp(-ret1)) < 1e-12);
+
double ret2 =
compute_state_update_lambert(h, unorm, L, s5_old[i][0]);
assert(std::abs(1.0 / h * ret2 - (s5_old[i] - unorm / L) / h -
exp(-ret2)) < 1e-12);
- assert(std::abs(ret2 - ret1) < 1e-14);
+
+ double ret3 =
+ compute_state_update_lambert_gsl(h, unorm, L, s4_old[i][0]);
+ assert(std::abs(1.0 / h * ret3 - (s4_old[i] - unorm / L) / h -
+ exp(-ret3)) < 1e-10);
+
+ assert(std::abs(ret1 - ret2) < 1e-14);
+ assert(std::abs(ret1 - ret3) < 1e-14);
+ assert(std::abs(ret2 - ret3) < 1e-14);
(*s5_new)[i][0] = ret1;
}
--
GitLab