From ad09881e5f2e16b5a5d33f373b2a5045d050ee12 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Thu, 2 Nov 2017 08:37:10 +0100 Subject: [PATCH] Add more Sacado::RFad and nested RFad-DFad examples to the testsuite These test the vector-mode functionality of the reverse AD numbers. --- tests/sacado/basic_01b.cc | 168 +++++++++++++++++++++++ tests/sacado/basic_01b.output | 8 ++ tests/sacado/basic_02b.cc | 251 ++++++++++++++++++++++++++++++++++ tests/sacado/basic_02b.output | 8 ++ 4 files changed, 435 insertions(+) create mode 100644 tests/sacado/basic_01b.cc create mode 100644 tests/sacado/basic_01b.output create mode 100644 tests/sacado/basic_02b.cc create mode 100644 tests/sacado/basic_02b.output diff --git a/tests/sacado/basic_01b.cc b/tests/sacado/basic_01b.cc new file mode 100644 index 0000000000..5f6c097a11 --- /dev/null +++ b/tests/sacado/basic_01b.cc @@ -0,0 +1,168 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Tests the computation of the first derivatives of a function using +// reverse mode AD with the Sacado::Rad::ADvar class. The difference between +// this test and sacado/basic_01a.cc is that we test the vector-mode capability +// of this number type. +// +// A related example that is shipped with Trilinos can be found at +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/tradvec_example.cpp + + +#include "../tests.h" + +#include + +// The functions to differentiate +template +NumberType f(const NumberType &x, const NumberType &y, const NumberType &z) +{ + return z*(x + z*y + x*y); +} +template +NumberType g(const NumberType &x, const NumberType &y, const NumberType &z) +{ + return std::sin(x*z)*std::cos(y/z); +} +template +NumberType h(const NumberType &x, const NumberType &y, const NumberType &z) +{ + return x*y*z; +} + +// The analytic derivative of the functions with respect to x and y +void +df(const double &x, const double &y, const double &z, + double &df_dx, double &df_dy) +{ + df_dx = z*(1.0 + y); + df_dy = z*(z + x); +} +void +dg(const double &x, const double &y, const double &z, + double &dg_dx, double &dg_dy) +{ + dg_dx = z*std::cos(x*z)*std::cos(y/z); + dg_dy = -(1.0/z)*std::sin(x*z)*std::sin(y/z); +} +void +dh(const double &x, const double &y, const double &z, + double &dh_dx, double &dh_dy) +{ + dh_dx = y*z; + dh_dy = x*z; +} + +int main() +{ + initlog(); + + // Values of function arguments + const double x = 5.0; + const double y = 10.0; + const double z = 4.0; + + // RAD objects: Independent variables + const Sacado::Rad::ADvar x_ad (x); + const Sacado::Rad::ADvar y_ad (y); + // RAD objects: Passive variables + const Sacado::Rad::ADvar z_ad (z); + + deallog << "x_ad: " << x_ad.val() << std::endl; + deallog << "y_ad: " << y_ad.val() << std::endl; + deallog << "z_ad: " << z_ad.val() << std::endl; + + // Compute functions + const double f = ::f(x, y, z); + const double g = ::g(x, y, z); + const double h = ::h(x, y, z); + + // Compute derivatives analytically + double df_dx = 0.0, df_dy = 0.0; + df(x, y, z, df_dx, df_dy); + double dg_dx = 0.0, dg_dy = 0.0; + dg(x, y, z, dg_dx, dg_dy); + double dh_dx = 0.0, dh_dy = 0.0; + dh(x, y, z, dh_dx, dh_dy); + + // Compute function values + // We specifically choose to do all of these computations + // before computing gradients, because this mixes the operations + // performed with each independent variables to produce each + // dependent variable + Sacado::Rad::ADvar f_rad = ::f(x_ad, y_ad, z_ad); // Cannot be const + Sacado::Rad::ADvar h_rad = ::h(x_ad, y_ad, z_ad); // Cannot be const <----- Before g_rad + Sacado::Rad::ADvar g_rad = ::g(x_ad, y_ad, z_ad); // Cannot be const + deallog << "f_rad: " << f_rad.val() << std::endl; + deallog << "g_rad: " << g_rad.val() << std::endl; + deallog << "h_rad: " << h_rad.val() << std::endl; + + // Configure the AD number to perform gradient computations + // related to the dependent function "f" + Sacado::Rad::ADvar::Outvar_Gradcomp(f_rad); + // Extract value and derivatives + const double f_ad = f_rad.val(); // f + const double df_dx_ad = x_ad.adj(); // df/dx + const double df_dy_ad = y_ad.adj(); // df/dy + + std::cout << "df_dx: " << df_dx << " df_dx_ad: " << df_dx_ad << std::endl; + std::cout << "df_dy: " << df_dy << " df_dy_ad: " << df_dy_ad << std::endl; + + // Configure the AD number to perform gradient computations + // related to the dependent function "g" + Sacado::Rad::ADvar::Outvar_Gradcomp(g_rad); + // Extract value and derivatives + const double g_ad = g_rad.val(); // g + const double dg_dx_ad = (x_ad.adj() - df_dx_ad); // dg/dx ; Note: Accumulation of partial derivatives + const double dg_dy_ad = (y_ad.adj() - df_dy_ad); // dg/dy ; Note: Accumulation of partial derivatives + + std::cout << "dg_dx: " << dg_dx << " dg_dx_ad: " << dg_dx_ad << std::endl; + std::cout << "dg_dy: " << dg_dy << " dg_dy_ad: " << dg_dy_ad << std::endl; + + // Configure the AD number to perform gradient computations + // related to the dependent function "h" + Sacado::Rad::ADvar::Outvar_Gradcomp(h_rad); + // Extract value and derivatives + const double h_ad = h_rad.val(); // h + const double dh_dx_ad = (x_ad.adj() - dg_dx_ad - df_dx_ad); // dh/dx ; Note: Accumulation of partial derivatives + const double dh_dy_ad = (y_ad.adj() - dg_dy_ad - df_dy_ad); // dh/dy ; Note: Accumulation of partial derivatives + // Observation: The accumulation of the adjoints appears to be related to + // the order in which ::Outvar_Gradcomp is called (i.e. which dependent + // variables the adjoints are computed for), rather than the order in + // which the functions themselves are evaluated. + + std::cout << "dh_dx: " << dh_dx << " dh_dx_ad: " << dh_dx_ad << std::endl; + std::cout << "dh_dy: " << dh_dy << " dh_dy_ad: " << dh_dy_ad << std::endl; + + const double tol = 1.0e-14; + Assert(std::fabs(f - f_ad) < tol, + ExcMessage("Computation incorrect: Value of f")); + Assert(std::fabs(df_dx - df_dx_ad) < tol && + std::fabs(df_dy - df_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative of f")); + Assert(std::fabs(g - g_ad) < tol, + ExcMessage("Computation incorrect: Value of g")); + Assert(std::fabs(dg_dx - dg_dx_ad) < tol && + std::fabs(dg_dy - dg_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative of g")); + Assert(std::fabs(h - h_ad) < tol, + ExcMessage("Computation incorrect: Value of h")); + Assert(std::fabs(dh_dx - dh_dx_ad) < tol && + std::fabs(dh_dy - dh_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative of h")); + + deallog << "OK" << std::endl; +} diff --git a/tests/sacado/basic_01b.output b/tests/sacado/basic_01b.output new file mode 100644 index 0000000000..28546c6220 --- /dev/null +++ b/tests/sacado/basic_01b.output @@ -0,0 +1,8 @@ + +DEAL::x_ad: 5.00000 +DEAL::y_ad: 10.0000 +DEAL::z_ad: 4.00000 +DEAL::f_rad: 380.000 +DEAL::g_rad: -0.731400 +DEAL::h_rad: 200.000 +DEAL::OK diff --git a/tests/sacado/basic_02b.cc b/tests/sacado/basic_02b.cc new file mode 100644 index 0000000000..9f167b267f --- /dev/null +++ b/tests/sacado/basic_02b.cc @@ -0,0 +1,251 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Tests the computation of the second derivatives of a function using +// nested reverse-forward mode AD. The Sacado::Rad::ADvar class is used to +// compute the first derivatives while the Sacado::Fad::DFad class, which utilizes +// dynamic memory allocation for the number of derivative components, is +// used for the second derivative calculations. The difference between +// this test and sacado/basic_01b.cc is that we test the vector-mode capability +// of Sacado::Rad::ADvar. +// +// A related examples that are shipped with Trilinos can be found at +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/trad_dfad_example.cpp +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/tradvec_example.cpp + + +#include "../tests.h" + +#include + +// The function to differentiate +template +NumberType +f(const NumberType &x, const NumberType &y, const NumberType2 &z) +{ + return z*(x*x*x + z*y*y + 0.5*x*y*y); +} +template +NumberType g(const NumberType &x, const NumberType &y, const NumberType &z) +{ + return std::sin(x*z)*std::cos(y/z); +} +template +NumberType h(const NumberType &x, const NumberType &y, const NumberType &z) +{ + return x*x*y*y*z; +} + +// The analytic derivative of the functions with respect to x and y +void +df(const double &x, const double &y, const double &z, + double &df_dx, double &df_dy) +{ + df_dx = z*(3.0*x*x + 0.5*y*y); + df_dy = z*(2.0*z*y + x*y); +} +void +dg(const double &x, const double &y, const double &z, + double &dg_dx, double &dg_dy) +{ + dg_dx = z*std::cos(x*z)*std::cos(y/z); + dg_dy = -(1.0/z)*std::sin(x*z)*std::sin(y/z); +} +void +dh(const double &x, const double &y, const double &z, + double &dh_dx, double &dh_dy) +{ + dh_dx = 2*x*y*y*z; + dh_dy = 2*x*x*y*z; +} + +// The analytic second derivatives of the functions with respect to x and y +void +d2f(const double &x, const double &y, const double &z, + double &d2f_dx_dx, double &d2f_dy_dy, + double &d2f_dy_dx) +{ + d2f_dx_dx = z*(6.0*x); + d2f_dy_dx = z*y; + d2f_dy_dy = z*(2.0*z + x); +} +void +d2g(const double &x, const double &y, const double &z, + double &d2g_dx_dx, double &d2g_dy_dy, + double &d2g_dy_dx) +{ + d2g_dx_dx = -z*z*std::sin(x*z)*std::cos(y/z); + d2g_dy_dx = -std::cos(x*z)*std::sin(y/z); + d2g_dy_dy = -(1.0/(z*z))*std::sin(x*z)*std::cos(y/z); +} +void +d2h(const double &x, const double &y, const double &z, + double &d2h_dx_dx, double &d2h_dy_dy, + double &d2h_dy_dx) +{ + d2h_dx_dx = 2*y*y*z; + d2h_dy_dx = 4*x*y*z; + d2h_dy_dy = 2*x*x*z; +} + +int main() +{ + initlog(); + + // Values of function arguments + const double x = -3.0; + const double y = 2.0; + const double z = 7.0; + + // Number of independent variables + const int num_deriv = 2; + + // FAD objects: Independent variables + Sacado::Rad::ADvar< Sacado::Fad::DFad > x_ad(Sacado::Fad::DFad(num_deriv, 0, x)); + Sacado::Rad::ADvar< Sacado::Fad::DFad > y_ad(Sacado::Fad::DFad(num_deriv, 1, y)); + // FAD objects: Passive variables + const Sacado::Rad::ADvar< Sacado::Fad::DFad > z_ad(z); + + deallog << "x_ad: " << x_ad.val() << std::endl; + deallog << "y_ad: " << y_ad.val() << std::endl; + deallog << "z_ad: " << z_ad.val() << std::endl; + + // Compute functions + const double f = ::f(x, y, z); + const double g = ::g(x, y, z); + const double h = ::h(x, y, z); + + // Compute derivatives analytically + double df_dx = 0.0, df_dy = 0.0; + df(x, y, z, df_dx, df_dy); + double dg_dx = 0.0, dg_dy = 0.0; + dg(x, y, z, dg_dx, dg_dy); + double dh_dx = 0.0, dh_dy = 0.0; + dh(x, y, z, dh_dx, dh_dy); + + // Compute second derivative analytically + double d2f_dx_dx = 0.0, d2f_dy_dy = 0.0, d2f_dy_dx = 0.0; + d2f(x, y, z, d2f_dx_dx, d2f_dy_dy, d2f_dy_dx); + double d2g_dx_dx = 0.0, d2g_dy_dy = 0.0, d2g_dy_dx = 0.0; + d2g(x, y, z, d2g_dx_dx, d2g_dy_dy, d2g_dy_dx); + double d2h_dx_dx = 0.0, d2h_dy_dy = 0.0, d2h_dy_dx = 0.0; + d2h(x, y, z, d2h_dx_dx, d2h_dy_dy, d2h_dy_dx); + + // Compute function values + // We specifically choose to do all of these computations + // before computing gradients, because this mixes the operations + // performed with each independent variables to produce each + // dependent variable + Sacado::Rad::ADvar< Sacado::Fad::DFad > f_rfad = ::f(x_ad, y_ad, z_ad); // Cannot be const + Sacado::Rad::ADvar< Sacado::Fad::DFad > h_rfad = ::h(x_ad, y_ad, z_ad); // Cannot be const <----- Before g_rad + Sacado::Rad::ADvar< Sacado::Fad::DFad > g_rfad = ::g(x_ad, y_ad, z_ad); // Cannot be const + deallog << "f_rfad: " << f_rfad.val() << std::endl; + deallog << "g_rfad: " << g_rfad.val() << std::endl; + deallog << "h_rfad: " << h_rfad.val() << std::endl; + + // Partial derivative accumulation terms + Sacado::Fad::DFad d_dx_rad_acc = 0.0; + Sacado::Fad::DFad d_dy_rad_acc = 0.0; + + // Configure the AD number to perform gradient computations + // related to the dependent function "f" + Sacado::Rad::ADvar< Sacado::Fad::DFad >::Outvar_Gradcomp(f_rfad); + // Extract value and derivatives + const double f_ad = f_rfad.val().val(); // f + const Sacado::Fad::DFad df_dx_fad = x_ad.adj(); // df/dx + const Sacado::Fad::DFad df_dy_fad = y_ad.adj(); // df/dy + const double df_dx_ad = df_dx_fad.val(); // df/dx + const double df_dy_ad = df_dy_fad.val(); // df/dy + const double d2f_dx_dx_ad = x_ad.adj().dx(0); // d^2f/dx^2 + const double d2f_dy_dx_ad = x_ad.adj().dx(1); // d^2f/dy_dx + const double d2f_dx_dy_ad = y_ad.adj().dx(0); // d^2f/dx_dy + const double d2f_dy_dy_ad = y_ad.adj().dx(1); // d^2f/dy^2 + + std::cout << "df_dx: " << df_dx << " df_dx_ad: " << df_dx_ad << std::endl; + std::cout << "df_dy: " << df_dy << " df_dy_ad: " << df_dy_ad << std::endl; + + // Configure the AD number to perform gradient computations + // related to the dependent function "g" + d_dx_rad_acc += df_dx_fad; + d_dy_rad_acc += df_dy_fad; + Sacado::Rad::ADvar< Sacado::Fad::DFad >::Outvar_Gradcomp(g_rfad); + // Extract value and derivatives + const double g_ad = g_rfad.val().val(); // g + const Sacado::Fad::DFad dg_dx_fad = x_ad.adj() - d_dx_rad_acc; // dg/dx ; Note: Accumulation of partial derivatives + const Sacado::Fad::DFad dg_dy_fad = y_ad.adj() - d_dy_rad_acc; // dg/dy ; Note: Accumulation of partial derivatives + const double dg_dx_ad = dg_dx_fad.val(); // dg/dx + const double dg_dy_ad = dg_dy_fad.val(); // dg/dy + const double d2g_dx_dx_ad = dg_dx_fad.dx(0); // d^2g/dx^2 + const double d2g_dy_dx_ad = dg_dx_fad.dx(1); // d^2g/dy_dx + const double d2g_dx_dy_ad = dg_dy_fad.dx(0); // d^2g/dx_dy + const double d2g_dy_dy_ad = dg_dy_fad.dx(1); // d^2g/dy^2 + + std::cout << "dg_dx: " << dg_dx << " dg_dx_ad: " << dg_dx_ad << std::endl; + std::cout << "dg_dy: " << dg_dy << " dg_dy_ad: " << dg_dy_ad << std::endl; + + // Configure the AD number to perform gradient computations + // related to the dependent function "h" + d_dx_rad_acc += dg_dx_fad; + d_dy_rad_acc += dg_dy_fad; + Sacado::Rad::ADvar< Sacado::Fad::DFad >::Outvar_Gradcomp(h_rfad); + // Extract value and derivatives + const double h_ad = h_rfad.val().val(); // h + const Sacado::Fad::DFad dh_dx_fad = x_ad.adj() - d_dx_rad_acc; // dh/dx ; Note: Accumulation of partial derivatives + const Sacado::Fad::DFad dh_dy_fad = y_ad.adj() - d_dy_rad_acc; // dh/dy ; Note: Accumulation of partial derivatives + const double dh_dx_ad = dh_dx_fad.val(); // dg/dx + const double dh_dy_ad = dh_dy_fad.val(); // dg/dy + const double d2h_dx_dx_ad = dh_dx_fad.dx(0); // d^2h/dx^2 + const double d2h_dy_dx_ad = dh_dx_fad.dx(1); // d^2h/dy_dx + const double d2h_dx_dy_ad = dh_dy_fad.dx(0); // d^2h/dx_dy + const double d2h_dy_dy_ad = dh_dy_fad.dx(1); // d^2h/dy^2 + // Observation: The accumulation of the adjoints appears to be related to + // the order in which ::Outvar_Gradcomp is called (i.e. which dependent + // variables the adjoints are computed for), rather than the order in + // which the functions themselves are evaluated. + + std::cout << "dh_dx: " << dh_dx << " dh_dx_ad: " << dh_dx_ad << std::endl; + std::cout << "dh_dy: " << dh_dy << " dh_dy_ad: " << dh_dy_ad << std::endl; + + const double tol = 1.0e-12; + Assert(std::fabs(f - f_ad) < tol, + ExcMessage("Computation incorrect: Value of f")); + Assert(std::fabs(df_dx - df_dx_ad) < tol && + std::fabs(df_dy - df_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative of f")); + Assert(std::fabs(d2f_dx_dx - d2f_dx_dx_ad) < tol && + std::fabs(d2f_dy_dy - d2f_dy_dy_ad) < tol && + std::fabs(d2f_dy_dx - d2f_dy_dx_ad) < tol, + ExcMessage("Computation incorrect: Second derivative of f")); + Assert(std::fabs(g - g_ad) < tol, + ExcMessage("Computation incorrect: Value of g")); + Assert(std::fabs(dg_dx - dg_dx_ad) < tol && + std::fabs(dg_dy - dg_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative of g")); + Assert(std::fabs(d2g_dx_dx - d2g_dx_dx_ad) < tol && + std::fabs(d2g_dy_dy - d2g_dy_dy_ad) < tol && + std::fabs(d2g_dy_dx - d2g_dy_dx_ad) < tol, + ExcMessage("Computation incorrect: Second derivative of g")); + Assert(std::fabs(h - h_ad) < tol, + ExcMessage("Computation incorrect: Value of h")); + Assert(std::fabs(dh_dx - dh_dx_ad) < tol && + std::fabs(dh_dy - dh_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative of h")); + Assert(std::fabs(d2h_dx_dx - d2h_dx_dx_ad) < tol && + std::fabs(d2h_dy_dy - d2h_dy_dy_ad) < tol && + std::fabs(d2h_dy_dx - d2h_dy_dx_ad) < tol, + ExcMessage("Computation incorrect: Second derivative of h")); + + deallog << "OK" << std::endl; +} diff --git a/tests/sacado/basic_02b.output b/tests/sacado/basic_02b.output new file mode 100644 index 0000000000..f8f4092dfc --- /dev/null +++ b/tests/sacado/basic_02b.output @@ -0,0 +1,8 @@ + +DEAL::x_ad: -3.00000 [ 1.00000 0.00000 ] +DEAL::y_ad: 2.00000 [ 0.00000 1.00000 ] +DEAL::z_ad: 7.00000 [ ] +DEAL::f_rfad: -35.0000 [ 203.000 154.000 ] +DEAL::g_rfad: -0.802738 [ -3.67867 0.0336865 ] +DEAL::h_rfad: 252.000 [ -168.000 252.000 ] +DEAL::OK -- 2.39.5