From d85186240e9f5a7e67255ff238b84ed3878b010f Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 6 Sep 2017 12:09:50 +0200 Subject: [PATCH] Add a quick test for Adol-C. --- tests/quick_tests/CMakeLists.txt | 5 ++ tests/quick_tests/adolc.cc | 111 +++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+) create mode 100644 tests/quick_tests/adolc.cc diff --git a/tests/quick_tests/CMakeLists.txt b/tests/quick_tests/CMakeLists.txt index e6a3529a02..678d2d6d07 100644 --- a/tests/quick_tests/CMakeLists.txt +++ b/tests/quick_tests/CMakeLists.txt @@ -146,6 +146,11 @@ IF (DEAL_II_WITH_ARPACK AND DEAL_II_WITH_UMFPACK) make_quicktest("arpack" ${_mybuild} "") ENDIF() +# Test Adol-C +IF (DEAL_II_WITH_ADOLC) + make_quicktest("adolc" ${_mybuild} "") +ENDIF() + # A custom test target: diff --git a/tests/quick_tests/adolc.cc b/tests/quick_tests/adolc.cc new file mode 100644 index 0000000000..bddd47e2f3 --- /dev/null +++ b/tests/quick_tests/adolc.cc @@ -0,0 +1,111 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test of basic functionality for Adol-C: +// - Taped doubles +// - Gradient and hessian computations +// Adapted from https://github.com/Homebrew/homebrew-science/blob/master/adol-c.rb + +#include + +#include +#include +#include + +#include + +using namespace dealii; + +int main(void) +{ + const double tol = 1e-12; + const unsigned int n = 10; + std::size_t tape_stats[STAT_SIZE]; + + double *xp = new double[n]; + double yp = 0.0; + adouble *x = new adouble[n]; + adouble y = 1.0; + + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + trace_on(1); + for (unsigned int i = 0; i < n; i++) + { + x[i] <<= xp[i]; + y *= x[i]; + } + y >>= yp; + delete[] x; + trace_off(); + tapestats(1, tape_stats); + + // --- Function --- + double *f = new double; + function(1, 1, n, xp, f); + + const double error_func_1 = yp - 1.0 / (1.0 + n); + const double error_func_2 = *f - 1.0 / (1.0 + n); + + deallog << "Error (function 1): " << error_func_1 << std::endl; + deallog << "Error (function 2): " << error_func_2 << std::endl; + Assert(error_func_1 < tol, ExcMessage("Should be zero!")); + Assert(error_func_2 < tol, ExcMessage("Should be zero!")); + + // --- Gradient --- + + double *g = new double[n]; + gradient(1, n, xp, g); + + double err_grad = 0; + for (unsigned int i = 0; i < n; i++) + err_grad += std::abs(g[i] - yp / xp[i]); + + deallog << "Error (gradient): " << err_grad << std::endl; + Assert(err_grad < tol, ExcMessage("Should be zero!")); + + // --- Hessian --- + + double **H = new double*[n]; + for (unsigned int i = 0; i < n; ++i) + H[i] = new double[i+1]; + + hessian(1, n, xp, H); + + double error_hess = 0; + for (unsigned int i = 0; i < n; i++) + for (unsigned int j = 0; j < n; j++) + if (i > j) + error_hess += std::abs(H[i][j] - g[i] / xp[j]); + + deallog << "Error (hessian): " << error_hess << std::endl; + Assert(error_hess < tol, ExcMessage("Should be zero!")); + + // -- Cleanup --- + + delete f; + f = nullptr; + + delete[] g; + g = nullptr; + + for (unsigned int i = 0; i < n; i++) + delete[] H[i]; + delete[] H; + H = nullptr; + + return 0; +} -- 2.39.5