From: Jean-Paul Pelteret Date: Fri, 11 Aug 2017 09:16:01 +0000 (-0600) Subject: Add stand-alone tests for Adol-C X-Git-Tag: v9.0.0-rc1~1278^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4783%2Fhead;p=dealii.git Add stand-alone tests for Adol-C --- diff --git a/tests/adolc/CMakeLists.txt b/tests/adolc/CMakeLists.txt new file mode 100644 index 0000000000..e3743f4f4e --- /dev/null +++ b/tests/adolc/CMakeLists.txt @@ -0,0 +1,6 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.9) +INCLUDE(../setup_testsubproject.cmake) +PROJECT(testsuite CXX) +IF(DEAL_II_WITH_ADOLC) + DEAL_II_PICKUP_TESTS() +ENDIF() diff --git a/tests/adolc/basic_function_gradient_hessian_taped.cc b/tests/adolc/basic_function_gradient_hessian_taped.cc new file mode 100644 index 0000000000..cd3d179eac --- /dev/null +++ b/tests/adolc/basic_function_gradient_hessian_taped.cc @@ -0,0 +1,105 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 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: +// - Taped doubles +// - Gradient and hessian computations +// Adapted from https://github.com/Homebrew/homebrew-science/blob/master/adol-c.rb + +#include "../tests.h" +#include +#include +#include + +#include + +int main(void) +{ + initlog(); + + 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; + + // --- 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; + + // --- 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; + + // -- 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; +} diff --git a/tests/adolc/basic_function_gradient_hessian_taped.output b/tests/adolc/basic_function_gradient_hessian_taped.output new file mode 100644 index 0000000000..1d6fd6ce4e --- /dev/null +++ b/tests/adolc/basic_function_gradient_hessian_taped.output @@ -0,0 +1,5 @@ + +DEAL::Error (function 1): 0 +DEAL::Error (function 2): 0 +DEAL::Error (gradient): 0 +DEAL::Error (hessian): 0 diff --git a/tests/adolc/basic_function_gradient_tapeless.cc b/tests/adolc/basic_function_gradient_tapeless.cc new file mode 100644 index 0000000000..73671902e1 --- /dev/null +++ b/tests/adolc/basic_function_gradient_tapeless.cc @@ -0,0 +1,65 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 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: +// - Tapeless doubles +// - Function and gradient computations +// Adapted from https://github.com/Homebrew/homebrew-science/blob/master/adol-c.rb +// See Adol-C manual v2.6 p65 figure 8 + +#include "../tests.h" +#include +#include + +#include + +int main(void) +{ + initlog(); + + const unsigned int n = 10; + adtl::setNumDir(n); + + adtl::adouble *x = new adtl::adouble[n]; + for (unsigned int i = 0; i < n; i++) + { + x[i] = (i + 1.0) / (2.0 + i); + x[i].setADValue(i,1); + } + + adtl::adouble y = 1.0; + for (unsigned int i = 0; i < n; i++) + y *= x[i]; + + // --- Function --- + + const double error_func = y.getValue() - 1.0 / (1.0 + n); + + deallog << "Error (function): " << error_func << std::endl; + + // --- Gradient --- + + double err_grad = 0; + for (unsigned int i = 0; i < n; i++) + err_grad += std::abs(y.getADValue(i) - y.getValue() / x[i].getValue()); + + deallog << "Error (gradient): " << err_grad << std::endl; + + // -- Cleanup --- + + delete[] x; + + return 0; +} diff --git a/tests/adolc/basic_function_gradient_tapeless.output b/tests/adolc/basic_function_gradient_tapeless.output new file mode 100644 index 0000000000..0218f0a198 --- /dev/null +++ b/tests/adolc/basic_function_gradient_tapeless.output @@ -0,0 +1,3 @@ + +DEAL::Error (function): 0 +DEAL::Error (gradient): 0 diff --git a/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc b/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc new file mode 100644 index 0000000000..9b3b2cb032 --- /dev/null +++ b/tests/adolc/basic_functions_jacobian_taped-cpp_vector.cc @@ -0,0 +1,117 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 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: +// - Taped doubles +// - Multiple dependent functions and Jacobian computations +// - Update to using C++ data structures where possible +// - Test that the tracing of active dependent/independent variables still +// works with these structures +// Adapted from https://github.com/Homebrew/homebrew-science/blob/master/adol-c.rb + +#include "../tests.h" +#include +#include +#include + +#include + +void test_reset_vector_values (const bool reset_values, const int tape_index) +{ + const unsigned int m = 5; // Dependents + const unsigned int n = 10; // Independents + + std::vector xp (n, 0.0); + std::vector yp (m, 0.0); + std::vector x (n, 1.0); // Dep. variable values initially set here + std::vector y (m, 1.0); + + if (reset_values == false) + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + trace_on(tape_index); + for (unsigned int i = 0; i < n; i++) + { + x[i] <<= xp[i]; + for (unsigned int j = 0; j < m; ++j) + y[j] *= (j+1)*x[i]; + } + for (unsigned int j = 0; j < m; ++j) + y[j] >>= yp[j]; + + trace_off(); + // tapestats(1, tape_stats); + + // --- Change values --- + if (reset_values == true) + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + // --- Functions --- + double *f = new double[m]; + function(tape_index, m, n, xp.data(), f); + + deallog << "Evaluation points:" << std::endl; + for (unsigned int i = 0; i < n; ++i) + deallog + << " x[" << i << "]: " << xp[i] + << std::endl; + + deallog << "Function values:" << std::endl; + for (unsigned int j = 0; j < m; ++j) + deallog + << " f[" << j << "]: " << f[j] + << " y[" << j << "]: " << yp[j] + << std::endl; + + // --- Jacobian --- + + double **J = new double*[m]; + for (unsigned int j = 0; j < m; ++j) + J[j] = new double[n]; + + jacobian(tape_index, m, n, xp.data(), J); + + deallog << "Function jacobian J:" << std::endl; + for (unsigned int j = 0; j < m; j++) + { + for (unsigned int i = 0; i < n; i++) + deallog << J[j][i] << (i +#include +#include + +#include + +int main(void) +{ + initlog(); + + const unsigned int m = 5; // Dependents + const unsigned int n = 10; // Independents + std::size_t tape_stats[STAT_SIZE]; + + double *xp = new double[n]; + double *yp = new double[m]; + adouble *x = new adouble[n]; + adouble *y = new adouble[m]; + + for (unsigned int i = 0; i < n; i++) + xp[i] = (i + 1.0) / (2.0 + i); + + for (unsigned int j = 0; j < m; ++j) + y[j] = 1.0; + + trace_on(1); + for (unsigned int i = 0; i < n; i++) + { + x[i] <<= xp[i]; + for (unsigned int j = 0; j < m; ++j) + y[j] *= (j+1)*x[i]; + } + for (unsigned int j = 0; j < m; ++j) + y[j] >>= yp[j]; + + delete[] x; + delete[] y; + + trace_off(); + tapestats(1, tape_stats); + + // --- Functions --- + double *f = new double[m]; + function(1, m, n, xp, f); + + deallog << "Evaluation points:" << std::endl; + for (unsigned int i = 0; i < n; ++i) + deallog + << " x[" << i << "]: " << xp[i] + << std::endl; + + deallog << "Function values:" << std::endl; + for (unsigned int j = 0; j < m; ++j) + deallog + << " f[" << j << "]: " << f[j] + << " y[" << j << "]: " << yp[j] + << std::endl; + + // --- Jacobian --- + + double **J = new double*[m]; + for (unsigned int j = 0; j < m; ++j) + J[j] = new double[n]; + + jacobian(1, m, n, xp, J); + + deallog << "Function jacobian J:" << std::endl; + for (unsigned int j = 0; j < m; j++) + { + for (unsigned int i = 0; i < n; i++) + deallog << J[j][i] << (i +#include + +#include + +int main(void) +{ + initlog(); + + const unsigned int m = 5; // Dependents + const unsigned int n = 10; // Independents + adtl::setNumDir(n); + + adtl::adouble *x = new adtl::adouble[n]; + for (unsigned int i = 0; i < n; i++) + { + x[i] = (i + 1.0) / (2.0 + i); + x[i].setADValue(i,1); + } + + adtl::adouble *y = new adtl::adouble[m]; + for (unsigned int j = 0; j < m; ++j) + y[j] = 1.0; + + for (unsigned int i = 0; i < n; i++) + for (unsigned int j = 0; j < m; ++j) + y[j] *= (j+1)*x[i]; + + // --- Functions --- + + deallog << "Evaluation points:" << std::endl; + for (unsigned int i = 0; i < n; ++i) + deallog + << " x[" << i << "]: " << x[i].getValue() + << std::endl; + + deallog << "Function values:" << std::endl; + for (unsigned int j = 0; j < m; ++j) + deallog + << " y[" << j << "]: " << y[j].getValue() + << std::endl; + + // --- Jacobian --- + + deallog << "Function jacobian J:" << std::endl; + for (unsigned int j = 0; j < m; j++) + { + for (unsigned int i = 0; i < n; i++) + deallog << y[j].getADValue(i) << (i