Merge branch 'master' into fadg

*AD integrator for pLaplacian
*Select between AD integrator and hond coded integrator

.gitignore

@@ -142,6 +142,16 @@ examples/petsc/mode_*
 
 examples/pumi/ex1
 examples/pumi/ex[126]p
+
+examples/hiop/ex9.mesh
+examples/hiop/ex9-mesh.*
+examples/hiop/ex9-init.*
+examples/hiop/ex9-final.*
+
+examples/ex71
+examples/ex71p
+examples/Example71*
+
 examples/pumi/refined.mesh
 examples/pumi/sol.gf
 examples/pumi/mesh.*

CMakeLists.txt

@@ -341,6 +341,34 @@ if (MFEM_USE_HIOP)
   # find_package updates HIOP_FOUND, HIOP_INCLUDE_DIRS, HIOP_LIBRARIES
 endif()
 
+# ADEPT package
+if (MFEM_USE_ADEPT)
+        find_package(ADEPT REQUIRED)
+        # find_package updates ADEPT_FOUND, ADEPT_INCLUDE_DIRS, ADEPT_LIBRARIES
+endif()
+
+# FADBAD++ package
+if (MFEM_USE_FADBADPP)
+	find_package(FADBADPP REQUIRED)
+        # find_package updates FADBADPP_FOUND, FADBADPP_INCLUDE_DIRS, FADBADPP_LIBRARIES
+endif()
+
+# CUDA
+if (MFEM_USE_CUDA)
+   set(CMAKE_CUDA_STANDARD 11)
+   set(CMAKE_CUDA_STANDARD_REQUIRED ON)
+   set(CMAKE_CUDA_EXTENSIONS OFF)
+   set(CMAKE_CUDA_FLAGS "-arch=${CUDA_ARCH} --expt-extended-lambda"
+       CACHE STRING "CUDA flags set for MFEM" FORCE)
+   if (MFEM_USE_MPI)
+      set(CUDA_CCBIN_COMPILER ${MPI_CXX_COMPILER})
+   else()
+      set(CUDA_CCBIN_COMPILER ${CMAKE_CXX_COMPILER})
+   endif()
+   string(APPEND CMAKE_CUDA_FLAGS " -ccbin ${CUDA_CCBIN_COMPILER}")
+   set(CMAKE_CUDA_HOST_LINK_LAUNCHER ${CUDA_CCBIN_COMPILER})
+endif()
+
 # OCCA
 if (MFEM_USE_OCCA)
    find_package(OCCA REQUIRED)
@@ -393,7 +421,8 @@ endif()
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     SLEPC MESQUITE SuperLUDist MUMPS STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
     MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2
-    CUSPARSE MKL_CPARDISO AMGX)
+    CUSPARSE MKL_CPARDISO AMGX FADBADPP)
+
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")

INSTALL

@@ -463,6 +463,19 @@ MFEM_USE_HIOP = YES/NO
   Enable the usage of HiOp (https://github.com/LLNL/hiop) in MFEM. HiOp is an
   HPC solver for nonlinear optimization problems.
 
+MFEM_USE_ADEPT = YES/NO
+  Enable automatic differentiation using the ADEPT library. 
+  (http://www.met.reading.ac.uk/clouds/adept) 
+  Please, compile the library with flag --disable-openmp.
+
+MFEM_USE_FADBADPP = YES/NO
+  Enable automatic differentiation using the FADBAD++ library.
+  www.fadbad.com/fadbad.html 
+
+MFEM_USE_ADFORWARD = YES/NO
+  Enable forward mode for AD packages. This option is valid
+  only if the AD package supports two modes (backward/forward).
+
 MFEM_USE_CUDA = YES/NO
    Enables support for CUDA devices in MFEM. CUDA is a parallel computing
    platform and programming model for general computing on graphical processing
@@ -674,6 +687,16 @@ The specific libraries and their options are:
   Options: HIOP_OPT, HIOP_LIB.
   Versions: HIOP >= 0.1.
 
+- ADEPT (optional), used with MFEM_USE_ADEPT = YES
+  URL: www.met.reading.ac.uk/clouds/adept/
+  Options: ADEPT_OPT, ADEPT_LIB
+  Versions: 1.1 and 2.0.5
+
+- FADBAD++ (optiobal), used with MFEM_USE_FADBADPP = YES
+  URL: www.fadbad.com/fadbad.html
+  Options: FADBADPP_OPT
+  Versions: 2.1
+
 - GSLIB (optional), used when MFEM_USE_GSLIB = YES. The gslib library must be
   built prior to the MFEM build, as follows: download gslib-1.0.5, untar it at
   the same level as MFEM and create a symbolic link: "ln -s gslib-1.0.5 gslib".
@@ -859,6 +882,9 @@ MFEM_USE_MPFR
 MFEM_USE_ZLIB
 MFEM_USE_PUMI
 MFEM_USE_HIOP
+MFEM_USE_ADEPT
+MFEM_USE_FADBADPP
+MFEM_USE_ADFORWARD
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_CEED
@@ -914,6 +940,8 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - POSIXCLOCKS
  - PUMI
  - HIOP
+ - ADEPT
+ - FADBAD++
  - OCCA
  - RAJA
  - UMPIRE

config/cmake/MFEMConfig.cmake.in

@@ -53,6 +53,9 @@ set(MFEM_USE_CEED @MFEM_USE_CEED@)
 set(MFEM_USE_UMPIRE @MFEM_USE_UMPIRE@)
 set(MFEM_USE_SIMD @MFEM_USE_SIMD@)
 set(MFEM_USE_ADIOS2 @MFEM_USE_ADIOS2@)
+set(MFEM_USE_ADEPT @MFEM_USE_ADEPT@)
+set(MFEM_USE_FADBADPP @MFEM_USE_FADBADPP@)
+set(MFEM_USE_ADFORWARD @MFEM_USE_ADFORWARD@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
 set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")

config/cmake/config.hpp.in

@@ -162,6 +162,16 @@
 // library.
 #cmakedefine MFEM_USE_SIMMETRIX
 
+
+// use ADEPT library for AD
+#cmakedefine MFEM_USE_ADEPT
+
+// use FADBAD++ library for AD
+#cmakedefine MFEM_USE_FADBADPP
+
+// use forward mode for automatic differentiation
+#cmakedefine MFEM_USE_ADFORWARD
+
 // Enable interface to the MKL CPardiso library.
 #cmakedefine MFEM_USE_MKL_CPARDISO
 

config/cmake/modules/FindADEPT.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Sets the following variables:
+#   - ADEPT_FOUND
+#   - ADEPT_INCLUDE_DIRS
+#   - ADEPT_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(ADEPT ADEPT ADEPT_DIR
+  "include" "adept.hpp"
+  "lib" "libadept.so"
+  "Paths to headers required by ADEPT."
+  "Libraries required by ADEPT.")
+

config/cmake/modules/FindFADBADPP.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Sets the following variables:
+#   - FADBADPP_FOUND
+#   - FADBADPP_INCLUDE_DIRS
+#   - FADBADPP_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(FADBADPP FADBADPP FADBADPP_DIR
+  "include" "fadiff.h"
+  "lib" ""
+  "Paths to headers required by FADBADPP."
+  "Libraries required by FADBADPP.")
+

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -742,7 +742,8 @@ function(mfem_export_mk_files)
       MFEM_USE_GNUTLS MFEM_USE_GSLIB MFEM_USE_NETCDF MFEM_USE_PETSC
       MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT MFEM_USE_PUMI
       MFEM_USE_CUDA MFEM_USE_OCCA MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD
-      MFEM_USE_ADIOS2)
+      MFEM_USE_ADIOS2 MFEM_USE_ADEPT MFEM_USE_FADBADPP MFEM_USE_ADFORWARD)
+
   foreach(var ${CONFIG_MK_BOOL_VARS})
     if (${var})
       set(${var} YES)

config/config.hpp.in

@@ -171,7 +171,18 @@
 // library.
 // #define MFEM_USE_SIMMETRIX
 
+
+// use ADEPT library for AD
+// #define MFEM_USE_ADEPT
+
+// use FADBAD++ library for AD
+// #define MFEM_USE_FADBADPP
+
+// use forward mode for automatic differentiation
+// #define MFEM_USE_ADFORWARD
+
 // Enable interface to the MKL CPardiso library.
 // #define MFEM_USE_MKL_CPARDISO
 
+
 #endif // MFEM_CONFIG_HEADER

config/config.mk.in

@@ -46,6 +46,9 @@ MFEM_USE_SIDRE         = @MFEM_USE_SIDRE@
 MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
 MFEM_USE_HIOP          = @MFEM_USE_HIOP@
+MFEM_USE_ADEPT         = @MFEM_USE_ADEPT@
+MFEM_USE_FADBADPP      = @MFEM_USE_FADBADPP@
+MFEM_USE_ADFORWARD     = @MFEM_USE_ADFORWARD@
 MFEM_USE_GSLIB         = @MFEM_USE_GSLIB@
 MFEM_USE_CUDA          = @MFEM_USE_CUDA@
 MFEM_USE_HIP           = @MFEM_USE_HIP@

config/defaults.cmake

@@ -55,8 +55,12 @@ option(MFEM_USE_CEED "Enable CEED" OFF)
 option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
 option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" OFF)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
+option(MFEM_USE_ADEPT "Enable AD using ADEPT" OFF)
+option(MFEM_USE_FADBADPP "Enable AD using FADBAD++" OFF)
+option(MFEM_USE_ADFORWARD "Enable forward mode for AD" OFF)
 option(MFEM_USE_MKL_CPARDISO "Enable MKL CPardiso" OFF)
 
+
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
 # Allow a user to disable testing, examples, and/or miniapps at CONFIGURE TIME
@@ -209,6 +213,13 @@ set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")
 set(LAPACK_LIBRARIES "" CACHE STRING "The LAPACK library.")
 
+set(ADEPT_INCLUDE_DIRS "${MFEM_DIR}/../adept-1.1/include" CACHE STRING "Path to ADEPT headers.")
+set(ADEPT_LIBRARIES "-L${MFEM_DIR}/../adept-1.1/lib -ladept" CACHE STRING "The ADEPT library.")
+
+set(FADBADPP_INCLUDE_DIRS "${MFEM_DIR}/../FADBAD++" CACHE STRING "Path to FADBAD++ headers.")
+set(FADBADPP_LIBRARIES "")
+
+
 # Some useful variables:
 set(CMAKE_SKIP_PREPROCESSED_SOURCE_RULES ON) # Skip *.i rules
 set(CMAKE_SKIP_ASSEMBLY_SOURCE_RULES  ON)    # Skip *.s rules

config/defaults.mk

@@ -145,6 +145,9 @@ MFEM_USE_CEED          = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
+MFEM_USE_ADEPT         = NO
+MFEM_USE_FADBADPP      = NO
+MFEM_USE_ADFORWARD     = NO
 MFEM_USE_MKL_CPARDISO  = NO
 
 # MPI library compile and link flags
@@ -374,6 +377,16 @@ HIOP_DIR = @MFEM_DIR@/../hiop/install
 HIOP_OPT = -I$(HIOP_DIR)/include
 HIOP_LIB = -L$(HIOP_DIR)/lib -lhiop $(LAPACK_LIB)
 
+# ADEPT
+ADEPT_DIR = @MFEM_DIR@/../adept-1.1
+ADEPT_OPT = -I$(ADEPT_DIR)/include
+ADEPT_LIB = -L$(ADEPT_DIR)/lib -ladept 
+
+# FADBAD++
+FADBADPP_DIR = @MFEM_DIR@/../FADBAD++
+FADBADPP_OPT = -I$(FADBADPP_DIR)
+FADBADPP_LIB = -L.
+
 # GSLIB library
 GSLIB_DIR = @MFEM_DIR@/../gslib/build
 GSLIB_OPT = -I$(GSLIB_DIR)/include

examples/CMakeLists.txt

@@ -34,6 +34,7 @@ list(APPEND ALL_EXE_SRCS
   ex25.cpp
   ex26.cpp
   ex27.cpp
+  ex71.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -64,6 +65,7 @@ if (MFEM_USE_MPI)
     ex25p.cpp
     ex26p.cpp
     ex27p.cpp
+    ex71p.cpp
     )
 endif()
 

examples/ex71.cpp

@@ -0,0 +1,456 @@
+//                       MFEM Example 71 - Serial Version
+//
+// Compile with: make ex71
+//
+// Sample runs:
+//    ex71 -m ../data/beam-quad.mesh -pp 3.5
+//    ex71 -m ../data/beam-tri.mesh  -pp 4.6
+//    ex71 -m ../data/beam-hex.mesh
+//    ex71 -m ../data/beam-tet.mesh
+//    ex71 -m ../data/beam-wedge.mesh
+//
+// Description:  This examples solves a quasi-static nonlinear
+//               p-Laplacian problem with zero Dirichlet boundary
+//               conditions applied on all defined boundaries
+//
+//           The example demonstrates the use of nonlinear operators combined
+//           with automatic differentiation (AD). The definitions of the
+//           integrators are written in the ex71.hpp. Selecting integrator=0
+//           will use the manually implemented integrator. Selecting
+//           integrator=1,2 will utilize one of the AD integrators.
+//
+//           The AD integrators are implemented in ex71.hpp (pLaplaceAD).
+//           The integrand qint is a function which is evaluated at every
+//           integration point. For implementations utilizing ADQFunctionTJ,
+//           the user has to implement the function and the residual
+//           evaluation. The Jacobian of the residual is evaluated using AD
+//
+//           For implementations utilizing ADQFunctionTH, the user has to
+//           implement only the function evaluation (as a template) and the
+//           first derivative (the residual) and the second derivatives (the
+//           Hessian) are evaluated using AD.
+//
+//           We recommend viewing examples 1 and 19, before viewing this
+//           example.
+
+#include "ex71.hpp"
+
+using namespace mfem;
+
+///Non-linear solver for the p-Laplacian problem.
+class NLSolverPLaplacian
+{
+public:
+   ///Constructor Input: imesh - FE mesh, finite element space,
+   /// power for the p-Laplacian, external load (source, input),
+   /// regularization parameter
+   NLSolverPLaplacian(Mesh& imesh, FiniteElementSpace& ifespace,
+                      double powerp=2,
+                      Coefficient* load=nullptr,
+                      double regularizationp=1e-7)
+   {
+      //default parameters for
+      //the Newton solver
+      newton_rtol = 1e-4;
+      newton_atol = 1e-8;
+      newton_iter = 10;
+
+      //linear solver
+      linear_rtol = 1e-7;
+      linear_atol = 1e-15;
+      linear_iter = 500;
+
+      print_level = 0;
+
+      //set the mesh
+      mesh=&imesh;
+
+      //set the fespace
+      fespace=&ifespace;
+
+      //set the parameters
+      plap_epsilon=new ConstantCoefficient(regularizationp);
+      plap_power=new ConstantCoefficient(powerp);
+      if (load==nullptr)
+      {
+         plap_input=new ConstantCoefficient(1.0);
+         input_ownership=true;
+      }
+      else
+      {
+         plap_input=load;
+         input_ownership=false;
+      }
+
+      //set the nonlinear form
+      nf=nullptr;
+      lsolver=nullptr;
+      prec=nullptr;
+      ns=nullptr;
+
+      //set the default integrator
+      integ=0; //hand coded
+
+   }
+
+   ~NLSolverPLaplacian()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) { delete ns;}
+      if (prec!=nullptr) { delete prec;}
+      if (lsolver!=nullptr) { delete lsolver;}
+      if (input_ownership) { delete plap_input;}
+      delete plap_epsilon;
+      delete plap_power;
+   }
+
+   ///Set the integrator.
+   /// 0 - hand coded, 1 - AD based (compute only Heassian by AD),
+   /// 2 - AD based (compute residual and Hessian by AD)
+   void SetIntegrator(int intr)
+   {
+      integ=intr;
+   }
+
+
+   //set relative tolerance for the Newton solver
+   void SetNRRTol(double rtol)
+   {
+      newton_rtol=rtol;
+   }
+
+   //set absolute tolerance for the Newton solver
+   void SetNRATol(double atol)
+   {
+      newton_atol=atol;
+   }
+
+   //set max iterations for the NR solver
+   void SetMaxNRIter(int miter)
+   {
+      newton_iter=miter;
+   }
+
+   void SetLSRTol(double rtol)
+   {
+      linear_rtol=rtol;
+   }
+
+   void SetLSATol(double atol)
+   {
+      linear_atol=atol;
+   }
+
+   //set max iterations for the linear solver
+   void SetMaxLSIter(int miter)
+   {
+      linear_iter=miter;
+   }
+
+   //set the print level
+   void SetPrintLevel(int plev)
+   {
+      print_level=plev;
+   }
+
+   ///The state vector is used as initial condition for the NR solver.
+   /// On return the statev holds the solution to the problem.
+   void Solve(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         AllocSolvers();
+      }
+      Vector b; //RHS is zero
+      ns->Mult(b, statev);
+   }
+
+   ///Compute the energy
+   double GetEnergy(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         //allocate the solvers
+         AllocSolvers();
+      }
+      return nf->GetEnergy(statev);
+   }
+
+
+private:
+
+   void AllocSolvers()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) {delete ns;}
+      if (prec!=nullptr) {delete prec;}
+      if (lsolver!=nullptr) { delete lsolver;}
+
+      // Define the essential boundary attributes
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      nf = new NonlinearForm(fespace);
+
+      if (integ==0)
+      {
+         nf->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else if (integ==1)
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<pLapIntegrandTJ>(*plap_power,
+                                                                 *plap_epsilon,*plap_input));
+      }
+      else
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<pLapIntegrandTH>(*plap_power,
+                                                                 *plap_epsilon,*plap_input));
+      }
+
+      nf->SetEssentialBC(ess_bdr);
+
+#ifdef MFEM_USE_SUITESPARSE
+      prec = new UMFPackSolver();
+#else
+      prec = new GSSmoother();
+#endif
+
+      //allocate the linear solver
+      lsolver=new CGSolver();
+      lsolver->SetRelTol(linear_rtol);
+      lsolver->SetAbsTol(linear_atol);
+      lsolver->SetMaxIter(linear_iter);
+      lsolver->SetPrintLevel(print_level);
+      lsolver->SetPreconditioner(*prec);
+
+      //allocate the NR solver
+      ns = new NewtonSolver();
+      ns->iterative_mode = true;
+      ns->SetSolver(*lsolver);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(newton_rtol);
+      ns->SetAbsTol(newton_atol);
+      ns->SetMaxIter(newton_iter);
+   }
+
+   double newton_rtol;
+   double newton_atol;
+   int newton_iter;
+
+   double linear_rtol;
+   double linear_atol;
+   int linear_iter;
+
+   int print_level;
+
+
+   //reference to the mesh
+   Mesh* mesh;
+   //reference to the fespace
+   FiniteElementSpace *fespace;
+
+   //nonlinear form for the p-laplacian
+   NonlinearForm *nf;
+   CGSolver *lsolver; //linear solver
+   Solver *prec; //preconditioner for the linear solver
+   NewtonSolver *ns; //NR solver
+   int integ;
+
+   //power of the p-laplacian
+   Coefficient* plap_power;
+   //regularization parammeter
+   Coefficient* plap_epsilon;
+   //load(input) paramater
+   Coefficient* plap_input;
+   bool input_ownership;
+};
+
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 3;
+   int order = 1;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 10;
+   int print_level = 0;
+
+   double pp = 2.0; // p-Lapalacian power
+
+   int integrator = 2; // 2 - use AD for Residual and Hessian
+   // 1 - use AD for Hessian only
+   // 0 - do not use AD (hand coded)
+   StopWatch *timer = new StopWatch();
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&newton_rel_tol,
+                  "-rel",
+                  "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol,
+                  "-abs",
+                  "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter,
+                  "-it",
+                  "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.AddOption(&pp,
+                  "-pp",
+                  "--power-parameter",
+                  "Power parameter (>=2.0) for the p-Laplacian.");
+   args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
+   args.AddOption(&integrator,
+                  "-int",
+                  "--integrator",
+                  "Integrator 0: standard; 1: AD for Hessian; 2: AD for residual and Hessian");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(std::cout);
+      return 1;
+   }
+   args.PrintOptions(std::cout);
+
+   // 2. Read the (serial) mesh from the given mesh file.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define the load parameter for the p-Laplacian
+   ConstantCoefficient load(1.00);
+
+   // 5. Define the finite element spaces for the solution
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(mesh, &fec, 1, Ordering::byVDIM);
+   int glob_size = fespace.GetTrueVSize();
+
+   std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+
+   // 6. Define the solution grid function
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 7. Define the solution true vector
+   Vector sv(fespace.GetTrueVSize());
+   sv = 0.0;
+
+   // 8. Define ParaView DataCollection
+   ParaViewDataCollection *dacol = new ParaViewDataCollection("Example71", mesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol", &x);
+
+   // 9. Define the NR solver
+   NLSolverPLaplacian* nr;
+
+   // 10. Start with linear diffusion - solvable for any initial guess
+   nr=new NLSolverPLaplacian(*mesh, fespace, 2.0, &load);
+   nr->SetIntegrator(integrator);
+   nr->SetMaxNRIter(newton_iter);
+   nr->SetNRATol(newton_abs_tol);
+   nr->SetNRRTol(newton_rel_tol);
+   timer->Clear();
+   timer->Start();
+   nr->Solve(sv);
+   timer->Stop();
+   std::cout << "[pp=2] The solution time is: " << timer->RealTime()
+             << std::endl;
+   // Compute the energy
+   double energy = nr->GetEnergy(sv);
+   std::cout << "[pp=2] The total energy of the system is E=" << energy
+             << std::endl;
+   delete nr;
+   x.SetFromTrueDofs(sv);
+   dacol->SetTime(2.0);
+   dacol->SetCycle(2);
+   dacol->Save();
+
+
+   // 11. Continue with powers higher than 2
+   for (int i = 3; i < pp; i++)
+   {
+      nr=new NLSolverPLaplacian(*mesh, fespace, (double)i, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(sv);
+      timer->Stop();
+      std::cout << "[pp=" << i
+                << "] The solution time is: " << timer->RealTime() << std::endl;
+      energy = nr->GetEnergy(sv);
+      std::cout << "[pp="<< i<<"] The total energy of the system is E=" << energy
+                << std::endl;
+      delete nr;
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(i);
+      dacol->SetCycle(i);
+      dacol->Save();
+   }
+
+   // 12. Continue with the final power
+   if (std::abs(pp - 2.0) > std::numeric_limits<double>::epsilon())
+   {
+      nr=new NLSolverPLaplacian(*mesh, fespace, pp, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(sv);
+      timer->Stop();
+      std::cout << "[pp=" << pp
+                << "] The solution time is: " << timer->RealTime() << std::endl;
+      energy = nr->GetEnergy(sv);
+      std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
+                << std::endl;
+      delete nr;
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(pp);
+      if (pp < 2.0)
+      {
+         dacol->SetCycle(std::floor(pp));
+      }
+      else
+      {
+         dacol->SetCycle(std::ceil(pp));
+      }
+      dacol->Save();
+   }
+
+   // 13. Free the memory
+   delete dacol;
+   delete mesh;
+   delete timer;
+   return 0;
+}

examples/ex71.hpp

@@ -0,0 +1,603 @@
+// Shared implementation ex71p/ex71 for the AD integrands and the manually
+// implemented integrators
+
+#ifndef EX71_HPP
+#define EX71_HPP
+
+#include "mfem.hpp"
+#include <memory>
+#include <iostream>
+#include <fstream>
+
+namespace mfem
+{
+
+///Example: Implementation of the energy and the residual
+/// for p-Laplacian problem. Both, the energy and the residual
+/// are evaluated at the integration points for PDE parameters
+/// vparam and state fields (derivatives with respect to x,y,z
+/// and primal field) stored in vector uu.
+template<typename DType, typename MVType>
+class MyQFunctorJ
+{
+public:
+   ///The operator returns the energy for the  p-Laplacian problem.
+   /// The input parameters vparam are: vparam[0] - the p-Laplacian
+   /// power, vparam[1] small value ensuring the exsitance of an unique
+   /// solution, and vparam[2] - the distributed extenal input to the
+   /// equation.
+   DType operator()(const Vector &vparam, MVType &uu)
+   {
+      double pp = vparam[0];
+      double ee = vparam[1];
+      double ff = vparam[2];
+
+      DType u = uu[3];
+      DType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
+
+      DType rez = pow(ee * ee + norm2, pp / 2.0) / pp - ff * u;
+      return rez;
+   }
+
+   ///The operator returns the first derivative of the energy with respect
+   /// to all state variables. These are set in vector uu and consist of the
+   /// derivatives with respect to x,y,z and the primal field. The derivative
+   /// is stored in vector rr with length equal to the length of vector uu.
+   void operator()(const Vector &vparam, MVType &uu, MVType &rr)
+   {
+      double pp = vparam[0];
+      double ee = vparam[1];
+      double ff = vparam[2];
+
+      DType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
+      DType tvar = pow(ee * ee + norm2, (pp - 2.0) / 2.0);
+
+      rr[0] = tvar * uu[0];
+      rr[1] = tvar * uu[1];
+      rr[2] = tvar * uu[2];
+      rr[3] = -ff;
+   }
+};
+
+///Defines AD class  pLapIntegrandTJ utilized for the automatic
+/// evaluation of the Hessian of the energy of the p-Laplacian.
+/// In general the length of the residual vector rr is not know
+/// in advance and therefore is provided explicitly as a second
+/// template argument in the definition of the class. The
+/// evaluation of the Hessian is based on the functor class MyQFunctorJ.
+typedef ADQFunctionTJ<MyQFunctorJ, 4> pLapIntegrandTJ;
+
+///Defines template class (functor) for evaluating the energy
+/// of the p-Laplacian problem. The input parameters vparam are:
+/// vparam[0] - the p-Laplacian power, vparam[1] small value
+/// ensuring exsitance of an unique solution, and vparam[2] -
+/// the distributed extenal input to the PDE.
+template<typename DType, typename MVType>
+class MyQFunctorH
+{
+public:
+   ///Returns the energy of a  p-Laplacian for state field input
+   /// provided in vector uu and parameters provided in vector
+   /// vparam.
+   DType operator()(const Vector &vparam, MVType &uu)
+   {
+      double pp = vparam[0];
+      double ee = vparam[1];
+      double ff = vparam[2];
+
+      DType u = uu[3];
+      DType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
+
+      DType rez = pow(ee * ee + norm2, pp / 2.0) / pp - ff * u;
+      return rez;
+   }
+};
+
+///Defines class pLapIntegrandTH for automatic evaluation
+/// of the first and second derivatives of the energy for
+/// a p-Laplacian problem. The energy is encoded in the
+/// operator()(...) of the template MyQFunctorH class.
+typedef ADQFunctionTH<MyQFunctorH> pLapIntegrandTH;
+
+//comment the line below in order to use
+//pLapIntegrandTJ for differentiation
+//the user interface for both TH and TJ versions
+//is exactly the same
+//#define USE_ADH
+
+
+///Implements integrator for a p-Laplacian problem.
+/// The integrator is based on a class QFunction utilized for
+/// evaluating the energy, the first derivative and the Hessian
+/// of the energy. QFunction shold be replaced with pLapIntegrandTH
+/// or pLapIntegrandTJ.
+template<class QFunction>
+class pLaplaceAD : public NonlinearFormIntegrator
+{
+protected:
+   Coefficient *pp;
+   Coefficient *coeff;
+   Coefficient *load;
+
+   QFunction qint;
+
+public:
+   pLaplaceAD()
+   {
+      coeff = nullptr;
+      pp = nullptr;
+   }
+
+   pLaplaceAD(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
+
+   pLaplaceAD(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
+      : pp(&pp_), coeff(&q), load(&ld_)
+   {}
+
+   virtual ~pLaplaceAD() {}
+
+   virtual double GetElementEnergy(const FiniteElement &el,
+                                   ElementTransformation &trans,
+                                   const Vector &elfun) override
+   {
+      double energy = 0.0;
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+
+      Vector vparam(3); //[power, epsilon, load]
+      Vector uu(4);     //[diff_x,diff_y,diff_z,u]
+
+      uu = 0.0;
+      vparam[0] = 2.0;  //default power
+      vparam[1] = 1e-8; //default epsilon
+      vparam[2] = 1.0;  //default load
+
+      double w;
+      double detJ;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         detJ = (square ? w : w * w);
+         w = ip.weight * w;
+
+         el.CalcDShape(ip, dshape_iso);
+         el.CalcShape(ip, shapef);
+         // AdjugateJacobian = / adj(J),         if J is square
+         //                    \ adj(J^t.J).J^t, otherwise
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided by detJ for obtaining the real value
+         // calculate the gradient
+         dshape_xyz.MultTranspose(elfun, grad);
+
+         // set the power
+         if (pp != nullptr)
+         {
+            vparam[0] = pp->Eval(trans, ip);
+         }
+
+         // set the coefficient ensuring positiveness of the tangent matrix
+         if (coeff != nullptr)
+         {
+            vparam[1] = coeff->Eval(trans, ip);
+         }
+         // add the contribution from the load
+         if (load != nullptr)
+         {
+            vparam[2] = load->Eval(trans, ip);
+         }
+         // fill the values of vector uu
+         for (int jj = 0; jj < spaceDim; jj++)
+         {
+            uu[jj] = grad[jj] / detJ;
+         }
+         uu[3] = shapef * elfun;
+
+         energy = energy + w * (qint.QFunction(vparam, uu));
+      }
+      return energy;
+   }
+
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &trans,
+                                      const Vector &elfun,
+                                      Vector &elvect) override
+   {
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector lvec(ndof);
+      elvect.SetSize(ndof);
+      elvect = 0.0;
+
+      DenseMatrix B(ndof, 4); //[diff_x,diff_y,diff_z, shape]
+      Vector vparam(3);       //[power, epsilon, load]
+      Vector uu(4);           //[diff_x,diff_y,diff_z,u]
+      Vector du(4);
+      B = 0.0;
+      uu = 0.0;
+      //initialize the parameters - keep the same order
+      //utilized in the pLapIntegrator definition
+      vparam[0] = 2.0;  //default power
+      vparam[1] = 1e-8; //default epsilon
+      vparam[2] = 1.0;  //default load
+
+      double w;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         //detJ = (square ? w : w*w);
+         w = ip.weight * w;
+
+         el.CalcDShape(ip, dshape_iso);
+         el.CalcShape(ip, shapef);
+         Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
+
+         // set the matrix B
+         for (int jj = 0; jj < spaceDim; jj++)
+         {
+            B.SetCol(jj, dshape_xyz.GetColumn(jj));
+         }
+         B.SetCol(3, shapef);
+
+         // set the power
+         if (pp != nullptr)
+         {
+            vparam[0] = pp->Eval(trans, ip);
+         }
+         // set the coefficient ensuring positiveness of the tangent matrix
+         if (coeff != nullptr)
+         {
+            vparam[1] = coeff->Eval(trans, ip);
+         }
+         // add the contribution from the load
+         if (load != nullptr)
+         {
+            vparam[2] = load->Eval(trans, ip);
+         }
+
+         // calculate uu
+         B.MultTranspose(elfun, uu);
+         // calculate derivative of the energy with respect to uu
+         qint.QFunctionDU(vparam, uu, du);
+
+         B.Mult(du, lvec);
+         elvect.Add(w, lvec);
+      } // end integration loop
+   }
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat) override
+   {
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      elmat.SetSize(ndof, ndof);
+      elmat = 0.0;
+
+      DenseMatrix B(ndof, 4); // [diff_x,diff_y,diff_z, shape]
+      DenseMatrix A(ndof, 4);
+      Vector vparam(3); // [power, epsilon, load]
+      Vector uu(4);     // [diff_x,diff_y,diff_z,u]
+      DenseMatrix duu(4, 4);
+      B = 0.0;
+      uu = 0.0;
+      // initialize the parameters - keep the same order
+      // utilized in the pLapIntegrator definition
+      vparam[0] = 2.0;  // default power
+      vparam[1] = 1e-8; // default epsilon
+      vparam[2] = 1.0;  // default load
+
+      double w;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         w = ip.weight * w;
+
+         el.CalcDShape(ip, dshape_iso);
+         el.CalcShape(ip, shapef);
+         Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
+
+         // set the matrix B
+         for (int jj = 0; jj < spaceDim; jj++)
+         {
+            B.SetCol(jj, dshape_xyz.GetColumn(jj));
+         }
+         B.SetCol(3, shapef);
+
+         // set the power
+         if (pp != nullptr)
+         {
+            vparam[0] = pp->Eval(trans, ip);
+         }
+         // set the coefficient ensuring positiveness of the tangent matrix
+         if (coeff != nullptr)
+         {
+            vparam[1] = coeff->Eval(trans, ip);
+         }
+         // add the contribution from the load
+         if (load != nullptr)
+         {
+            vparam[2] = load->Eval(trans, ip);
+         }
+
+         // calculate uu
+         B.MultTranspose(elfun, uu);
+         // calculate derivative of the energy with respect to uu
+         qint.QFunctionDD(vparam, uu, duu);
+
+         Mult(B, duu, A);
+         AddMult_a_ABt(w, A, B, elmat);
+
+      } // end integration loop
+   }
+};
+
+///Implements hand-coded integrator for a p-Laplacian problem.
+/// Utilized as alternative for the  pLaplaceAD class based on
+/// automatic differentiation.
+class pLaplace : public NonlinearFormIntegrator
+{
+protected:
+   Coefficient *pp;
+   Coefficient *coeff;
+   Coefficient *load;
+
+public:
+   pLaplace()
+   {
+      coeff = nullptr;
+      pp = nullptr;
+   }
+
+   pLaplace(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
+
+   pLaplace(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
+      : pp(&pp_), coeff(&q), load(&ld_)
+   {}
+
+   virtual ~pLaplace() {}
+
+   virtual double GetElementEnergy(const FiniteElement &el,
+                                   ElementTransformation &trans,
+                                   const Vector &elfun) override
+   {
+      double energy = 0.0;
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+
+      double w;
+      double detJ;
+      double nrgrad2;
+      double ppp = 2.0;
+      double eee = 0.0;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         detJ = (square ? w : w * w);
+         w = ip.weight * w;
+
+         el.CalcDShape(ip, dshape_iso);
+         el.CalcShape(ip, shapef);
+         // AdjugateJacobian = / adj(J),         if J is square
+         //                    \ adj(J^t.J).J^t, otherwise
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided by detJ for obtaining the real value
+         // calculate the gradient
+         dshape_xyz.MultTranspose(elfun, grad);
+         nrgrad2 = grad * grad / (detJ * detJ);
+
+         // set the power
+         if (pp != nullptr)
+         {
+            ppp = pp->Eval(trans, ip);
+         }
+
+         // set the coefficient ensuring positiveness of the tangent matrix
+         if (coeff != nullptr)
+         {
+            eee = coeff->Eval(trans, ip);
+         }
+
+         energy = energy + w * std::pow(nrgrad2 + eee * eee, ppp / 2.0) / ppp;
+
+         // add the contribution from the load
+         if (load != nullptr)
+         {
+            energy = energy - w * (shapef * elfun) * load->Eval(trans, ip);
+         }
+      }
+      return energy;
+   }
+
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &trans,
+                                      const Vector &elfun,
+                                      Vector &elvect) override
+   {
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+      Vector lvec(ndof);
+      elvect.SetSize(ndof);
+      elvect = 0.0;
+
+      double w;
+      double detJ;
+      double nrgrad;
+      double aa;
+      double ppp = 2.0;
+      double eee = 0.0;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         detJ = (square ? w : w * w);
+         w = ip.weight * w; //w;
+
+         el.CalcDShape(ip, dshape_iso);
+         el.CalcShape(ip, shapef);
+         // AdjugateJacobian = / adj(J),         if J is square
+         //                    \ adj(J^t.J).J^t, otherwise
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided by detJ for obtaining the real value
+
+         // calculate the gradient
+         dshape_xyz.MultTranspose(elfun, grad);
+         nrgrad = grad.Norml2() / detJ;
+         // grad is not scaled so far, i.e., grad=grad/detJ
+
+         // set the power
+         if (pp != nullptr)
+         {
+            ppp = pp->Eval(trans, ip);
+         }
+
+         // set the coefficient ensuring positiveness of the tangent matrix
+         if (coeff != nullptr)
+         {
+            eee = coeff->Eval(trans, ip);
+         }
+
+         aa = nrgrad * nrgrad + eee * eee;
+         aa = std::pow(aa, (ppp - 2.0) / 2.0);
+         dshape_xyz.Mult(grad, lvec);
+         elvect.Add(w * aa / (detJ * detJ), lvec);
+
+         // add loading
+         if (load != nullptr)
+         {
+            elvect.Add(-w * load->Eval(trans, ip), shapef);
+         }
+      } // end integration loop
+   }
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat) override
+   {
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+      Vector lvec(ndof);
+      elmat.SetSize(ndof, ndof);
+      elmat = 0.0;
+
+      double w;
+      double detJ;
+      double nrgrad;
+      double aa0;
+      double aa1;
+      double ppp = 2.0;
+      double eee = 0.0;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         detJ = (square ? w : w * w);
+         w = ip.weight * w;
+
+         el.CalcDShape(ip, dshape_iso);
+         // AdjugateJacobian = / adj(J),         if J is square
+         //                    \ adj(J^t.J).J^t, otherwise
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided by detJ for obtaining the real value
+         // grad is not scaled so far,i.e., grad=grad/detJ
+
+         //set the power
+         if (pp != nullptr)
+         {
+            ppp = pp->Eval(trans, ip);
+         }
+         //set the coefficient ensuring positiveness of the tangent matrix
+         if (coeff != nullptr)
+         {
+            eee = coeff->Eval(trans, ip);
+         }
+
+         //calculate the gradient
+         dshape_xyz.MultTranspose(elfun, grad);
+         nrgrad = grad.Norml2() / detJ;
+         aa0 = nrgrad * nrgrad + eee * eee;
+         aa1 = std::pow(aa0, (ppp - 2.0) / 2.0);
+         aa0 = (ppp - 2.0) * std::pow(aa0, (ppp - 4.0) / 2.0);
+         dshape_xyz.Mult(grad, lvec);
+         w = w / (detJ * detJ);
+         AddMult_a_VVt(w * aa0 / (detJ * detJ), lvec, elmat);
+         AddMult_a_AAt(w * aa1, dshape_xyz, elmat);
+
+      } // end integration loop
+   }
+};
+
+} // namespace mfem
+#endif

examples/ex71p.cpp

@@ -0,0 +1,505 @@
+//                       MFEM Example 71 - Parallel Version
+//
+// Compile with: make ex71p
+//
+// Sample runs:
+//   mpirun -np 2 ex71p -m ../data/beam-quad.mesh -pp 3.8
+//   mpirun -np 2 ex71p -m ../data/beam-tri.mesh  -pp 7.2
+//   mpirun -np 2 ex71p -m ../data/beam-hex.mesh
+//   mpirun -np 2 ex71p -m ../data/beam-tet.mesh
+//   mpirun -np 2 ex71p -m ../data/beam-wedge.mesh
+//
+// Description:  This examples solves a quasi-static nonlinear
+//               p-Laplacian problem with zero Dirichlet boundary
+//               conditions applied on all defined boundaries
+//
+//           The example demonstrates the use of nonlinear operators combined
+//           with automatic differentiation (AD). The definitions of the
+//           integrators are written in the ex71.hpp.  Selecting integrator=0
+//           will use the manually implemented integrator.  Selecting
+//           integrator=1,2 will utilize the AD integrator.
+//
+//        The AD integrators are implemented in ex71.hpp (pLaplaceAD).
+//           The integrand qint is a function which is evaluated at every
+//           integration point. For implementations utilizing ADQFunctionTJ,
+//           the user has to implement the function and the residual
+//           evaluation. The Jacobian of the residual is evaluated using AD
+//           For implementations utilizing ADQFunctionTH, the user has to
+//           implement only the function evaluation (as a template) and the
+//           first derivative (the residual) and the second derivatives (the
+//           Hessian) are evaluated using AD.
+//
+//           We recommend viewing examples 1 and 19, before viewing this
+//           example.
+
+#include "ex71.hpp"
+
+using namespace mfem;
+
+///Non-linear solver for the p-Laplacian problem.
+class ParNLSolverPLaplacian
+{
+public:
+   ///Constructor Input: imesh - FE mesh, finite element space,
+   /// power for the p-Laplacian, external load (source, input),
+   /// regularization parameter
+   ParNLSolverPLaplacian(MPI_Comm comm, ParMesh& imesh,
+                         ParFiniteElementSpace& ifespace,
+                         double powerp=2,
+                         Coefficient* load=nullptr,
+                         double regularizationp=1e-7)
+   {
+      lcomm = comm;
+
+      //default parameters for
+      //the Newton solver
+      newton_rtol = 1e-4;
+      newton_atol = 1e-8;
+      newton_iter = 10;
+
+      //linear solver
+      linear_rtol = 1e-7;
+      linear_atol = 1e-15;
+      linear_iter = 500;
+
+      print_level = 0;
+
+      //set the mesh
+      mesh=&imesh;
+
+      //set the fespace
+      fespace=&ifespace;
+
+      //set the parameters
+      plap_epsilon=new ConstantCoefficient(regularizationp);
+      plap_power=new ConstantCoefficient(powerp);
+      if (load==nullptr)
+      {
+         plap_input=new ConstantCoefficient(1.0);
+         input_ownership=true;
+      }
+      else
+      {
+         plap_input=load;
+         input_ownership=false;
+      }
+
+      nf=nullptr;
+      ns=nullptr;
+      gmres=nullptr;
+      prec=nullptr;
+
+      //set the default integrator
+      integ=0; //hand coded
+   }
+
+   ~ParNLSolverPLaplacian()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) { delete ns;}
+      if (prec!=nullptr) { delete prec;}
+      if (gmres!=nullptr) { delete gmres;}
+      if (input_ownership) { delete plap_input;}
+      delete plap_epsilon;
+      delete plap_power;
+   }
+
+   ///Set the integrator.
+   /// 0 - hand coded, 1 - AD based (compute only Heassian by AD),
+   /// 2 - AD based (compute residual and Hessian by AD)
+   void SetIntegrator(int intr)
+   {
+      integ=intr;
+   }
+
+   //set relative tolerance for the Newton solver
+   void SetNRRTol(double rtol)
+   {
+      newton_rtol=rtol;
+   }
+
+   //set absolute tolerance for the Newton solver
+   void SetNRATol(double atol)
+   {
+      newton_atol=atol;
+   }
+
+   //set max iterations for the NR solver
+   void SetMaxNRIter(int miter)
+   {
+      newton_iter=miter;
+   }
+
+   void SetLSRTol(double rtol)
+   {
+      linear_rtol=rtol;
+   }
+
+   void SetLSATol(double atol)
+   {
+      linear_atol=atol;
+   }
+
+   //set max iterations for the linear solver
+   void SetMaxLSIter(int miter)
+   {
+      linear_iter=miter;
+   }
+
+   //set the print level
+   void SetPrintLevel(int plev)
+   {
+      print_level=plev;
+   }
+
+   ///The state vector is used as initial condition for the NR solver.
+   /// On return the statev holds the solution to the problem.
+   void Solve(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         AllocSolvers();
+      }
+      Vector b; //RHS is zero
+      ns->Mult(b, statev);
+   }
+
+   ///Compute the energy
+   double GetEnergy(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         //allocate the solvers
+         AllocSolvers();
+      }
+      return nf->GetEnergy(statev);
+   }
+
+private:
+   void AllocSolvers()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) { delete ns;}
+      if (gmres!=nullptr) { delete gmres;}
+      if (prec!=nullptr) { delete prec;}
+
+      // Define the essential boundary attributes
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      nf = new ParNonlinearForm(fespace);
+      if (integ==0)
+      {
+         nf->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else if (integ==1)
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<pLapIntegrandTJ>(*plap_power,
+                                                                 *plap_epsilon,*plap_input));
+      }
+      else
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<pLapIntegrandTH>(*plap_power,
+                                                                 *plap_epsilon,*plap_input));
+      }
+
+      nf->SetEssentialBC(ess_bdr);
+
+      prec = new HypreBoomerAMG();
+      prec->SetPrintLevel(print_level);
+
+      gmres = new GMRESSolver(lcomm);
+      gmres->SetAbsTol(linear_atol);
+      gmres->SetRelTol(linear_rtol);
+      gmres->SetMaxIter(linear_iter);
+      gmres->SetPrintLevel(print_level);
+      gmres->SetPreconditioner(*prec);
+
+      ns = new NewtonSolver(lcomm);
+
+      ns->iterative_mode = true;
+      ns->SetSolver(*gmres);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(newton_rtol);
+      ns->SetAbsTol(newton_atol);
+      ns->SetMaxIter(newton_iter);
+   }
+
+   double newton_rtol;
+   double newton_atol;
+   int newton_iter;
+
+   double linear_rtol;
+   double linear_atol;
+   int linear_iter;
+
+   int print_level;
+
+   //power of the p-laplacian
+   Coefficient* plap_power;
+   //regularization parammeter
+   Coefficient* plap_epsilon;
+   //load(input) paramater
+   Coefficient* plap_input;
+   bool input_ownership;
+
+   MPI_Comm lcomm;
+
+   ParMesh *mesh;
+   ParFiniteElementSpace *fespace;
+
+   ParNonlinearForm *nf;
+
+   HypreBoomerAMG *prec;
+   GMRESSolver *gmres;
+   NewtonSolver *ns;
+   int integ;
+
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI
+   int num_procs, myrank;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+   // 2. Parse command-line options
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 3;
+   int par_ref_levels = 1;
+   int order = 2;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 10;
+   int print_level = 0;
+   double pp = 2.0;
+   int integrator = 2; //use AD
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels,
+                  "-rp",
+                  "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&newton_rel_tol,
+                  "-rel",
+                  "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol,
+                  "-abs",
+                  "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter,
+                  "-it",
+                  "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.AddOption(&pp,
+                  "-pp",
+                  "--power-parameter",
+                  "Power parameter (>=2.0) for the p-Laplacian.");
+   args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
+   args.AddOption(&integrator,
+                  "-int",
+                  "--integrator",
+                  "Integrator 0: standard; 1: AD for Hessian; 2: AD for residual and Hessian");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(std::cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   StopWatch *timer = new StopWatch();
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
+   //    with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define the load for the p-Laplacian
+   ConstantCoefficient load(1.00);
+
+   // 7. Define the finite element spaces for the solution
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(pmesh, &fec, 1, Ordering::byVDIM);
+   HYPRE_Int glob_size = fespace.GlobalTrueVSize();
+   if (myrank == 0)
+   {
+      std::cout << "Number of finite element unknowns: " << glob_size
+                << std::endl;
+   }
+
+   // 8. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+   HypreParVector *sv = x.GetTrueDofs();
+
+   // 9. Define ParaView DataCollection
+   ParaViewDataCollection *dacol = new ParaViewDataCollection("Example71",
+                                                              pmesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol", &x);
+
+   // 10. Define the NR solver
+   ParNLSolverPLaplacian* nr;
+
+   // 11. Start with linear diffusion - solvable for any initial guess
+   nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, 2.0, &load);
+   nr->SetIntegrator(integrator);
+   nr->SetMaxNRIter(newton_iter);
+   nr->SetNRATol(newton_abs_tol);
+   nr->SetNRRTol(newton_rel_tol);
+   nr->SetPrintLevel(print_level);
+   timer->Clear();
+   timer->Start();
+   nr->Solve(*sv);
+   timer->Stop();
+   if (myrank==0)
+   {
+      std::cout << "[pp=2] The solution time is: " << timer->RealTime()
+                << std::endl;
+   }
+   // Compute the energy
+   double energy = nr->GetEnergy(*sv);
+   if (myrank==0)
+   {
+      std::cout << "[pp=2] The total energy of the system is E=" << energy
+                << std::endl;
+   }
+   delete nr;
+   x.SetFromTrueDofs(*sv);
+   dacol->SetTime(2.0);
+   dacol->SetCycle(2);
+   dacol->Save();
+
+
+   // 12. Continue with powers higher than 2
+   for (int i = 3; i < pp; i++)
+   {
+      nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, (double)i, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      nr->SetPrintLevel(print_level);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<i<<"] The solution time is: " << timer->RealTime()
+                   << std::endl;
+      }
+      // Compute the energy
+      double energy = nr->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<i<<"] The total energy of the system is E=" << energy
+                   << std::endl;
+      }
+      delete nr;
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime((double)i);
+      dacol->SetCycle(i);
+      dacol->Save();
+   }
+
+   // 13. Continue with the final power
+   if (std::abs(pp - 2.0) > std::numeric_limits<double>::epsilon())
+   {
+      nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, pp, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      nr->SetPrintLevel(print_level);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<pp<<"] The solution time is: " << timer->RealTime()
+                   << std::endl;
+      }
+      // Compute the energy
+      double energy = nr->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
+                   << std::endl;
+      }
+      delete nr;
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime(pp);
+      if (pp < 2.0)
+      {
+         dacol->SetCycle(std::floor(pp));
+      }
+      else
+      {
+         dacol->SetCycle(std::ceil(pp));
+      }
+      dacol->Save();
+   }
+
+   // 14. Free the used memory
+   delete dacol;
+   delete sv;
+   delete pmesh;
+   delete timer;
+
+   MPI_Finalize();
+
+   return 0;
+}

examples/makefile

@@ -22,10 +22,10 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27
+  ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex71
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
   ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p ex25p\
-  ex26p ex27p
+  ex26p ex27p ex71p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -157,3 +157,4 @@ clean-exec:
 	@rm -f ex21*.mesh ex21*.sol ex21p_*.*
 	@rm -f ex23.mesh ex23-*.gf
 	@rm -f ex25.mesh ex25-*.gf ex25p-*.*
+	@rm -rf Example71

fem/CMakeLists.txt

@@ -110,6 +110,7 @@ set(HDRS
   tmop.hpp
   tmop_tools.hpp
   gslib.hpp
+  adnonlininteg.hpp
   transfer.hpp
   )
 

fem/adnonlininteg.hpp

@@ -0,0 +1,507 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_ADNONLININTEG
+#define MFEM_ADNONLININTEG
+
+#include "../config/config.hpp"
+#include "fe.hpp"
+#include "coefficient.hpp"
+#include "fespace.hpp"
+#include "nonlininteg.hpp"
+#include "../linalg/tadvector.hpp"
+#include "../linalg/taddensemat.hpp"
+#include "../linalg/fdual.hpp"
+
+#if defined MFEM_USE_ADEPT
+#include <adept.h>
+#elif defined MFEM_USE_FADBADPP
+#include <fadiff.h>
+#include <badiff.h>
+#endif
+
+
+namespace mfem
+{
+
+///Automatic differentiation class - for user coded functor CTD returning  PDE's
+/// energy and residual at a point, provides the derivative of the residual
+/// with respect to the active arguments of the functor.
+///
+/**
+   The two templated parameters are
+   template<typename, typename> class CTD and int m.
+   The templated class CTD is a user-provided functor which implements
+   the function of interest. The CTD class should have the following signature:
+
+template<typename DType, typename MVType>
+
+class MyQFunctorJ
+
+{
+
+    DType operator()(const Vector &vparam, MVType &uu)
+
+        {...}
+
+    void operator()(const Vector &vparam, MVType &uu, MVType &rr)
+
+    {...}
+
+};
+
+
+DType - represents the scalars in the functor, and  MVType - represents
+the vectors in the functor. In ADQFunctionTJ,  DType will be replaced
+either with double or AD-type.
+The first operator()(const Vector &vparam, MVType &uu) takes as input
+standard MFEM Vector vparam (passive arguments), which holds all parameters
+ supplied to the functor and MVType vector uu, which holds all active
+arguments provided to the functor. The operator returns a scalar
+representing the value of the function(energy) for passive parameters vparam
+and active arguments uu.
+
+
+The second void operator()(const Vector &vparam, MVType &uu, MVType &rr)
+returns a vector-valued function rr with passive input arguments vparam
+and active (AD) arguments supplied in vector uu. Since the size of the
+return vector rr is not known in advance, it should be provided explicitly
+to the ADQFunctionTJ class, i.e., the template parameter m in ADQFunctionTJ
+represents the dimensions of the return vector rr.
+
+
+For PDE discretized problem, the first operator should return the energy
+evaluated a point, and the second operator should return the residual
+evaluated at the point. The derivative of the residual vector rr with
+respect to the active arguments uu is evaluated automatically by
+ADQFunctionTJ and returned by the ADQFunctionTJ:: QFunctionDD(..) method.
+*/
+template<template<typename, typename> class CTD, int m>
+class ADQFunctionTJ
+{
+   // m - dimension of the residual vector
+   // the Jacobian will have dimensions [m,length(uu)]
+protected:
+#ifdef MFEM_USE_ADEPT
+   adept::Stack m_stack;
+#endif
+
+public:
+#if defined MFEM_USE_ADEPT
+   /// AD-type based on the Adept AD library.
+   typedef adept::adouble ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#elif defined MFEM_USE_FADBADPP
+#ifdef MFEM_USE_ADFORWARD
+   /// AD-type based on the FADBAD++ library -
+   /// used only for forward differentiation.
+   typedef fadbad::F<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#else
+   /// AD-type based on the FADBAD++ library -
+   /// used only for reverse AD mode.
+   typedef fadbad::B<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#endif
+#else
+   /// MFEM native forward AD-type
+   typedef ad::FDual<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#endif
+
+#ifdef MFEM_USE_ADEPT
+   ADQFunctionTJ() : m_stack(false) {}
+#else
+   ADQFunctionTJ() {}
+#endif
+
+   ~ADQFunctionTJ() {}
+
+   ///Returns the energy for passive arguments vparam and
+   /// active arguments uu. The evaluation is based on the
+   /// first operator in the user-supplied CTD template class.
+   double QFunction(const Vector &vparam, Vector &uu)
+   {
+      CTD<double, Vector> func;
+      return func(vparam, uu);
+   }
+
+   ///Returns vector valued function rr for passive input parametrs vparam
+   /// and active arguments provided in uu.
+   void QFunctionDU(const Vector &vparam, ADFVector &uu, ADFVector &rr)
+   {
+      CTD<ADFType, ADFVector> func;
+      func(vparam, uu, rr);
+   }
+
+   ///Evaluates automatically the first derivative of
+   /// QFunction(const Vector &vparam, Vector &uu) with respect to
+   /// the active vector arguments uu. The dimension or rr is the same
+   /// as the dimension of uu. The method can be used for testing
+   /// the correctness of hand-coded gradients of  QFunction(...).
+   void QFunctionAU(const Vector &vparam, Vector &uu, Vector &rr)
+   {
+      //the result is computed automatically by differentiating
+      //QFunction with respect to uu
+      CTD<ADFType, ADFVector> func;
+      int n = uu.Size();
+      rr.SetSize(n);
+
+#if defined MFEM_USE_ADEPT
+      //use ADEPT package
+      adept::Stack *p_stack = adept::active_stack();
+      p_stack->deactivate();
+
+      m_stack.activate();
+      {
+         ADFVector aduu(uu);
+         ADFType rez;
+         m_stack.new_recording();
+         rez = func(vparam, aduu);
+         m_stack.independent(aduu.GetData(), n); //independent variables
+         m_stack.dependent(&rez, 1);             //dependent variables
+         m_stack.jacobian(rr.GetData());
+      }
+      m_stack.deactivate();
+#elif defined MFEM_USE_FADBADPP
+      //use FADBAD++
+#ifdef MFEM_USE_ADFORWARD
+      {
+         ADFVector aduu(uu);
+         ADFType rez;
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].diff(ii, n);
+         }
+         rez = func(vparam, aduu);
+         for (int ii = 0; ii < n; ii++)
+         {
+            rr[ii] = rez.d(ii);
+         }
+      }
+#else
+      {
+         ADFVector aduu(uu);
+         ADFType rez;
+         rez = func(vparam, aduu);
+         rez.diff(0, 1);
+         for (int ii = 0; ii < n; ii++)
+         {
+            rr[ii] = aduu[ii].d(0);
+         }
+      }
+#endif
+#else
+      //use native AD package
+      {
+         ADFVector aduu(uu); //all dual numbers are initialized to zero
+         ADFType rez;
+
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].dual(1.0);
+            rez = func(vparam, aduu);
+            rr[ii] = rez.dual();
+            aduu[ii].dual(0.0);
+         }
+      }
+#endif
+   }
+
+   ///Returns vector valued function rr for supplied passive arguments
+   /// vparam and active arguments uu. The evaluation is based on the
+   /// user supplied CTD template class (second operator).
+   void QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+   {
+      CTD<double, Vector> func;
+      func(vparam, uu, rr);
+   }
+
+   ///Evaluates automatically the derivative or the vector function
+   /// QFunctionDU(...), i.e., the retuned vector rr, with respect to
+   /// the active arguments uu. The dimensions of the the dense matix
+   /// jac are [m,n] where m is the size of the vector rr and n is the
+   /// size of the active vector uu. The parameter m should be supplied as
+   /// template parameter int m in ADQFunctionTJ.
+   void QFunctionDD(const Vector &vparam, Vector &uu, DenseMatrix &jac)
+   {
+#if defined MFEM_USE_ADEPT
+      //use ADEPT package
+      adept::Stack *p_stack = adept::active_stack();
+      p_stack->deactivate();
+
+      int n = uu.Size();
+      jac.SetSize(m, n);
+      jac = 0.0;
+      m_stack.activate();
+      {
+         ADFVector aduu(uu);
+         ADFVector rr(m); //residual vector
+         m_stack.new_recording();
+         QFunctionDU(vparam, aduu, rr);
+         m_stack.independent(aduu.GetData(), n); //independent variables
+         m_stack.dependent(rr.GetData(), m);     //dependent variables
+         m_stack.jacobian(jac.Data());
+      }
+      m_stack.deactivate();
+#elif defined MFEM_USE_FADBADPP
+      //use FADBAD++
+#ifdef MFEM_USE_ADFORWARD
+      int n = uu.Size();
+      jac.SetSize(m, n);
+      jac = 0.0;
+      {
+         ADFVector aduu(uu);
+         ADFVector rr(m);
+
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].diff(ii, n);
+         }
+         QFunctionDU(vparam, aduu, rr);
+         for (int ii = 0; ii < n; ii++)
+         {
+            for (int jj = 0; jj < m; jj++)
+            {
+               jac(jj, ii) = rr[jj].d(ii);
+            }
+         }
+      }
+#else
+      int n = uu.Size();
+      jac.SetSize(m, n);
+      jac = 0.0;
+      {
+         ADFVector aduu(uu);
+         ADFVector rr(m);
+         QFunctionDU(vparam, aduu, rr);
+         for (int ii = 0; ii < m; ii++)
+         {
+            rr[ii].diff(ii, m);
+         }
+         for (int ii = 0; ii < n; ii++)
+         {
+            for (int jj = 0; jj < m; jj++)
+            {
+               jac(jj, ii) = aduu[ii].d(jj);
+            }
+         }
+      }
+#endif
+#else
+      //use native AD package
+      int n = uu.Size();
+      jac.SetSize(m, n);
+      jac = 0.0;
+      {
+         ADFVector aduu(uu); //all dual numbers are initialized to zero
+         ADFVector rr(m);
+
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].dual(1.0);
+            QFunctionDU(vparam, aduu, rr);
+            for (int jj = 0; jj < m; jj++)
+            {
+               jac(jj, ii) = rr[jj].dual();
+            }
+            aduu[ii].dual(0.0);
+         }
+      }
+#endif
+   }
+};
+
+//template class for differentiation; the function
+//for differentiation is supplied as a functor
+//the operator()(scalar,vector) defines the actual function
+
+///ADQFunctionTH is a templated class evaluating the first and
+/// the second derivatives of a user supplied function encoded
+/// in a functor class CTD. The signature of CTD is as follows:
+/**
+
+template<typename DType, typename MVType> class CTD
+
+{
+
+    DType operator()(const Vector &vparam, MVType &uu)
+    {...}
+
+};
+
+The operator returns a scalar representing the value of the
+function for passive arguments vparam and active arguments uu.
+The first derivative of the operator is evaluated automatically
+by ADQFunctionTH::QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+method. The length of the return vector rr is the same as for
+the input vector uu.  The second derivate (the Hessian) of the function
+is evaluated again automatically by
+QFunctionDD(const Vector &vparam, const Vector &uu, DenseMatrix &jac).
+
+Hessian evaluation is an expensive process. The provided AD-functionality
+is intended to be utilized for development purposes. The performance will
+increase significantly by replacing the AD-provided derivatives with hand-coded.
+*/
+template<template<typename, typename> class CTD>
+class ADQFunctionTH
+{
+public:
+#if defined MFEM_USE_FADBADPP
+   ///AD-type derived from FADBAD++
+   typedef fadbad::B<double> ADFType;
+   ///AD vector type for evaluation of first derivatives
+   typedef TADVector<ADFType> ADFVector;
+   ///AD dense matrix type for evaluation of first derivatives
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+   ///AD-type for the second derivatives derived from
+   ///FADBAD++
+   typedef fadbad::B<fadbad::F<double>> ADSType;
+   /// AD vector type for evaluation of second derivatives
+   typedef TADVector<ADSType> ADSVector;
+   /// AD dense matrix type for  evaluation of second derivatives
+   typedef TADDenseMatrix<ADSType> ADSDenseMatrix;
+#else
+   ///MFEM native AD-type for first derivatives
+   typedef ad::FDual<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+   ///MFEM native AD-type for second derivatives
+   typedef ad::FDual<ADFType> ADSType;
+   typedef TADVector<ADSType> ADSVector;
+   typedef TADDenseMatrix<ADSType> ADSDenseMatrix;
+#endif
+
+   ADQFunctionTH() {}
+
+   ~ADQFunctionTH() {}
+
+   ///Evaluates a function for arguments vparam and uu.
+   /// The evaluatin is based on the operator() in the
+   /// user provided functor CTD.
+   double QFunction(const Vector &vparam, Vector &uu)
+   {
+      CTD<double, Vector> tf;
+      return tf(vparam, uu);
+   }
+
+   ///Evaluates the first derivative of QFunction(...).
+   /// Intended for internal use only.
+   ADFType QFunction(const Vector &vparam, ADFVector &uu)
+   {
+      CTD<ADFType, ADFVector> tf;
+      return tf(vparam, uu);
+   }
+
+   ///Evaluates the second derivative of QFunction(...).
+   /// Intended for internal use only.
+   ADSType QFunction(const Vector &vparam, ADSVector &uu)
+   {
+      CTD<ADSType, ADSVector> tf;
+      return tf(vparam, uu);
+   }
+
+   ///Returns the first derivative of QFunction(...) with
+   /// respect to the active arguments proved in vector uu.
+   /// The length of rr is the same as for uu.
+   void QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+   {
+#if defined MFEM_USE_FADBADPP
+      int n = uu.Size();
+      rr.SetSize(n);
+      ADFVector aduu(uu);
+      ADFType rez;
+      rez = QFunction(vparam, aduu);
+      rez.diff(0, 1);
+      for (int ii = 0; ii < n; ii++)
+      {
+         rr[ii] = aduu[ii].d(0);
+      }
+#else
+      int n = uu.Size();
+      rr.SetSize(n);
+      ADFVector aduu(uu);
+      ADFType rez;
+      for (int ii = 0; ii < n; ii++)
+      {
+         aduu[ii].dual(1.0);
+         rez = QFunction(vparam, aduu);
+         rr[ii] = rez.dual();
+         aduu[ii].dual(0.0);
+      }
+#endif
+   }
+
+   ///Returns the Hessian of QFunction(...) in the dense matrix jac.
+   /// The dimensions of jac are m x m, where m is the length of vector uu.
+   void QFunctionDD(const Vector &vparam, const Vector &uu, DenseMatrix &jac)
+   {
+#if defined MFEM_USE_FADBADPP
+      int n = uu.Size();
+      jac.SetSize(n);
+      jac = 0.0;
+      {
+         ADSVector aduu(n);
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii] = uu[ii];
+            aduu[ii].x().diff(ii, n);
+         }
+         ADSType rez = QFunction(vparam, aduu);
+         rez.diff(0, 1);
+         for (int ii = 0; ii < n; ii++)
+         {
+            for (int jj = 0; jj < ii; jj++)
+            {
+               jac(ii, jj) = aduu[ii].d(0).d(jj);
+               jac(jj, ii) = aduu[jj].d(0).d(ii);
+            }
+            jac(ii, ii) = aduu[ii].d(0).d(ii);
+         }
+      }
+#else
+      int n = uu.Size();
+      jac.SetSize(n);
+      jac = 0.0;
+      {
+         ADSVector aduu(n);
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].real(ADFType(uu[ii], 0.0));
+            aduu[ii].dual(ADFType(0.0, 0.0));
+         }
+
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].real(ADFType(uu[ii], 1.0));
+            for (int jj = 0; jj < (ii + 1); jj++)
+            {
+               aduu[jj].dual(ADFType(1.0, 0.0));
+               ADSType rez = QFunction(vparam, aduu);
+               jac(ii, jj) = rez.dual().dual();
+               jac(jj, ii) = rez.dual().dual();
+               aduu[jj].dual(ADFType(0.0, 0.0));
+            }
+            aduu[ii].real(ADFType(uu[ii], 0.0));
+         }
+      }
+#endif
+   }
+
+};
+
+} // namespace mfem
+
+#endif

fem/fem.hpp

@@ -35,6 +35,7 @@
 #include "tmop.hpp"
 #include "tmop_tools.hpp"
 #include "gslib.hpp"
+#include "adnonlininteg.hpp"
 #include "restriction.hpp"
 #include "quadinterpolator.hpp"
 #include "quadinterpolator_face.hpp"

fem/nonlinearform.cpp

@@ -578,7 +578,9 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
          }
       }
 
+
    // free the allocated memory
+
    for (int i = 0; i < fes.Size(); ++i)
    {
       delete el_x[i];

linalg/fdual.hpp

@@ -0,0 +1,548 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef FDUAL_H
+#define FDUAL_H
+
+#include <cmath>
+#include <type_traits>
+
+namespace mfem
+{
+namespace ad
+{
+/** The FDual template class provides forward automatic differentiation
+   implementation based on dual numbers. The derivative of an arbitrary
+ function double f(double a)  can be obtained by replacing the double
+ type for the return value and the argument a  with FDual<double>, i.e.,
+ FDual<double> f(FDual<double> a). The derivative is evaluated automatically
+ by calling the function r=f(a). The value of the function is stored in r.pr
+ and the derivative in r.du. These can be extracted by the corresponding
+  methods real()/prim() and dual(). Internally, the function f can be
+ composed of standard functions predefined for FDual type. These consist
+ of a large set of functions replicating the functionality of the standard
+ math library, i.e., sin, cos, exp, log, ...
+
+New functions (non-member) can be eaily added to the class. Example:
+
+template<typename tbase>
+inline FDual<tbase> cos(const FDual<tbase> &f)
+{
+   return FDual<tbase>(cos(f.real()), -f.dual() * sin(f.real()));
+}
+
+The real part of the return value consists of the standard real value
+ of the function, i.e., cos(f.real()).
+
+The dual part of the return value consists of the first derivative of
+ the function with respect to the real part of the argument -sin(f.reaf)
+ multiplied with the dual part of the argument f.dual().
+*/
+template<typename tbase>
+class FDual
+{
+private:
+   /// Real value
+   tbase pr;
+   /// Dual value holding derivative information
+   tbase du;
+
+public:
+   /// Standard contructor - both values are set to zero.
+   FDual() : pr(0), du(0) {}
+
+   /// The constructor utilized in nested definition of dual numbers.
+   /// It is used for second and higher order derivatives.
+   template<class fltyp,
+            class = typename std::enable_if<std::is_arithmetic<fltyp>::value>::type>
+   FDual(fltyp &f) : pr(f), du(0)
+   {}
+
+   template<class fltyp,
+            class = typename std::enable_if<std::is_arithmetic<fltyp>::value>::type>
+   FDual(const fltyp &f) : pr(f), du(0)
+   {}
+
+   /// Standard constructor with user supplied inpur for both parts of the dual number.
+   FDual(tbase &pr_, tbase &du_) : pr(pr_), du(du_) {}
+
+   FDual(const tbase &pr_, const tbase &du_) : pr(pr_), du(du_) {}
+
+   FDual(FDual<tbase> &nm) : pr(nm.pr), du(nm.du) {}
+
+   FDual(const FDual<tbase> &nm) : pr(nm.pr), du(nm.du) {}
+
+   /// Return the real value of the dual number.
+   tbase prim() const { return pr; }
+
+   /// Same as prim(). Return the real value of the dual number.
+   tbase real() const { return pr; }
+
+   /// Return the dual value of the dual number.
+   tbase dual() const { return du; }
+
+   /// Set the primal and the dual values.
+   void set(const tbase &pr_, const tbase &du_)
+   {
+      pr = pr_;
+      du = du_;
+   }
+
+   void prim(const tbase &pr_) { pr = pr_; }
+
+   void real(const tbase &pr_) { pr = pr_; }
+
+   void dual(const tbase &du_) { du = du_; }
+
+   FDual<tbase> &operator=(tbase sc_)
+   {
+      pr = sc_;
+      du = tbase(0);
+      return *this;
+   }
+
+   FDual<tbase> &operator+=(tbase sc_)
+   {
+      pr = pr + sc_;
+      return *this;
+   }
+
+   FDual<tbase> &operator-=(tbase sc_)
+   {
+      pr = pr - sc_;
+      return *this;
+   }
+
+   FDual<tbase> &operator*=(tbase sc_)
+   {
+      pr = pr * sc_;
+      du = du * sc_;
+      return *this;
+   }
+
+   FDual<tbase> &operator/=(tbase sc_)
+   {
+      pr = pr / sc_;
+      du = du / sc_;
+      return *this;
+   }
+
+   FDual<tbase> &operator=(const FDual<tbase> &f)
+   {
+      pr = f.real();
+      du = f.dual();
+      return *this;
+   }
+
+   FDual<tbase> &operator+=(const FDual<tbase> &f)
+   {
+      pr += f.real();
+      du += f.dual();
+      return *this;
+   }
+
+   FDual<tbase> &operator-=(const FDual<tbase> &f)
+   {
+      pr -= f.real();
+      du -= f.dual();
+      return *this;
+   }
+
+   FDual<tbase> &operator*=(const FDual<tbase> &f)
+   {
+      du = du * f.real();
+      du = du + pr * f.dual();
+      pr = pr * f.real();
+      return *this;
+   }
+
+   FDual<tbase> &operator/=(const FDual<tbase> &f_)
+   {
+      pr = pr / f_.real();
+      du = du - pr * f_.dual();
+      du = du / f_.real();
+      return *this;
+   }
+};
+
+// non-member functions
+// boolean operations
+template<typename tbase>
+inline bool operator==(const FDual<tbase> &a1, const FDual<tbase> &a2)
+{
+   return a1.real() == a2.real();
+}
+
+template<typename tbase>
+inline bool operator==(tbase a, const FDual<tbase> &f_)
+{
+   return a == f_.real();
+}
+
+template<typename tbase>
+inline bool operator==(const FDual<tbase> &a, tbase b)
+{
+   return a.real() == b;
+}
+
+template<typename tbase>
+inline bool operator<(const FDual<tbase> &f1, const FDual<tbase> &f2)
+{
+   return f1.real() < f2.real();
+}
+
+template<typename tbase>
+inline bool operator<(const FDual<tbase> &f, tbase a)
+{
+   return f.real() < a;
+}
+
+template<typename tbase>
+inline bool operator<(tbase a, const FDual<tbase> &f)
+{
+   return a < f.real();
+}
+
+template<typename tbase>
+inline bool operator>(const FDual<tbase> &f1, const FDual<tbase> &f2)
+{
+   return f1.real() > f2.real();
+}
+
+template<typename tbase>
+inline bool operator>(const FDual<tbase> &f, tbase a)
+{
+   return f.real() > a;
+}
+
+template<typename tbase>
+inline bool operator>(tbase a, const FDual<tbase> &f)
+{
+   return (a > f.real());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator-(const FDual<tbase> &f)
+{
+   return FDual<tbase>(-f.real(), -f.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator-(const FDual<tbase> &f, tbase a)
+{
+   return FDual<tbase>(f.real() - a, f.dual());
+}
+
+template<typename tbase>
+inline FDual<FDual<tbase>> operator-(const FDual<FDual<tbase>> &f, tbase a)
+{
+   return FDual<FDual<tbase>>(f.real() - a, f.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator+(const FDual<tbase> &f, tbase a)
+{
+   return FDual<tbase>(f.real() + a, f.dual());
+}
+
+template<typename tbase>
+inline FDual<FDual<tbase>> operator+(const FDual<FDual<tbase>> &f, tbase a)
+{
+   return FDual<FDual<tbase>>(f.real() + a, f.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator*(const FDual<tbase> &f, tbase a)
+{
+   return FDual<tbase>(f.real() * a, f.dual() * a);
+}
+
+template<typename tbase>
+inline FDual<tbase> operator/(const FDual<tbase> &f, tbase a)
+{
+   return FDual<tbase>(f.real() / a, f.dual() / a);
+}
+
+template<typename tbase>
+inline FDual<FDual<tbase>> operator/(const FDual<FDual<tbase>> &f, tbase a)
+{
+   return FDual<FDual<tbase>>(f.real() / a, f.dual() / a);
+}
+
+template<typename tbase>
+inline FDual<tbase> operator+(tbase a, const FDual<tbase> &f)
+{
+   return FDual<tbase>(a + f.real(), f.dual());
+}
+
+template<typename tbase>
+inline FDual<FDual<tbase>> operator+(tbase a, const FDual<FDual<tbase>> &f)
+{
+   return FDual<FDual<tbase>>(a + f.real(), f.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator-(tbase a, const FDual<tbase> &f)
+{
+   return FDual<tbase>(a - f.real(), -f.dual());
+}
+
+template<typename tbase>
+inline FDual<FDual<tbase>> operator-(tbase a, const FDual<FDual<tbase>> &f)
+{
+   return FDual<FDual<tbase>>(a - f.real(), -f.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator*(tbase a, const FDual<tbase> &f)
+{
+   return FDual<tbase>(f.real() * a, f.dual() * a);
+}
+
+template<typename tbase>
+inline FDual<FDual<tbase>> operator*(tbase a, const FDual<FDual<tbase>> &f)
+{
+   return FDual<FDual<tbase>>(f.real() * a, f.dual() * a);
+}
+
+template<typename tbase>
+inline FDual<tbase> operator/(tbase a, const FDual<tbase> &f)
+{
+   a = a / f.real();
+   return FDual<tbase>(a, -a * f.dual() / f.real());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator+(const FDual<tbase> &f1, const FDual<tbase> &f2)
+{
+   return FDual<tbase>(f1.real() + f2.real(), f1.dual() + f2.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator-(const FDual<tbase> &f1, const FDual<tbase> &f2)
+{
+   return FDual<tbase>(f1.real() - f2.real(), f1.dual() - f2.dual());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator*(const FDual<tbase> &f1, const FDual<tbase> &f2)
+{
+   return FDual<tbase>(f1.real() * f2.real(),
+                       f1.real() * f2.dual() + f1.dual() * f2.real());
+}
+
+template<typename tbase>
+inline FDual<tbase> operator/(const FDual<tbase> &f1, const FDual<tbase> &f2)
+{
+   tbase a = tbase(1) / f2.real();
+   tbase b = f1.real() * a;
+   return FDual<tbase>(b, (f1.dual() - f2.dual() * b) * a);
+}
+
+template<typename tbase>
+inline FDual<tbase> acos(const FDual<tbase> &f)
+{
+   return FDual<tbase>(acos(f.real()),
+                       -f.dual() / sqrt(tbase(1) - f.real() * f.real()));
+}
+
+template<>
+inline FDual<double> acos(const FDual<double> &f)
+{
+   return FDual<double>(std::acos(f.real()),
+                        -f.dual() / std::sqrt(double(1) - f.real() * f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> asin(const FDual<tbase> &f)
+{
+   return FDual<tbase>(asin(f.real()),
+                       f.dual() / sqrt(tbase(1) - f.real() * f.real()));
+}
+
+template<>
+inline FDual<double> asin(const FDual<double> &f)
+{
+   return FDual<double>(std::asin(f.real()),
+                        f.dual() / std::sqrt(double(1) - f.real() * f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> atan(const FDual<tbase> &f)
+{
+   return FDual<tbase>(atan(f.real()),
+                       f.dual() / (tbase(1) + f.real() * f.real()));
+}
+
+template<>
+inline FDual<double> atan(const FDual<double> &f)
+{
+   return FDual<double>(std::atan(f.real()),
+                        f.dual() / (double(1) + f.real() * f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> cos(const FDual<tbase> &f)
+{
+   return FDual<tbase>(cos(f.real()), -f.dual() * sin(f.real()));
+}
+
+template<>
+inline FDual<double> cos(const FDual<double> &f)
+{
+   return FDual<double>(std::cos(f.real()), -f.dual() * std::sin(f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> cosh(const FDual<tbase> &f)
+{
+   return FDual<tbase>(cosh(f.real()), f.dual() * sinh(f.real()));
+}
+
+template<>
+inline FDual<double> cosh(const FDual<double> &f)
+{
+   return FDual<double>(std::cosh(f.real()), f.dual() * std::sinh(f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> exp(const FDual<tbase> &f)
+{
+   tbase x = exp(f.real());
+   return FDual<tbase>(x, f.dual() * x);
+}
+
+template<>
+inline FDual<double> exp(const FDual<double> &f)
+{
+   double x = std::exp(f.real());
+   return FDual<double>(x, f.dual() * x);
+}
+
+template<typename tbase>
+inline FDual<tbase> log(const FDual<tbase> &f)
+{
+   return FDual<tbase>(log(f.real()), f.dual() / f.real());
+}
+
+template<>
+inline FDual<double> log(const FDual<double> &f)
+{
+   return FDual<double>(std::log(f.real()), f.dual() / f.real());
+}
+
+template<typename tbase>
+inline FDual<tbase> log10(const FDual<tbase> &f)
+{
+   return log(f) / log(tbase(10));
+}
+
+template<>
+inline FDual<double> log10(const FDual<double> &f)
+{
+   return log(f) / std::log(double(10));
+}
+
+template<typename tbase>
+inline FDual<tbase> pow(const FDual<tbase> &a, const FDual<tbase> &b)
+{
+   return exp(log(a) * b);
+}
+
+template<typename tbase, typename tbase1>
+inline FDual<tbase> pow(const FDual<tbase> &a, const tbase1 &b)
+{
+   return exp(log(a) * tbase(b));
+}
+
+template<typename tbase, typename tbase1>
+inline FDual<tbase> pow(const tbase1 &a, const FDual<tbase> &b)
+{
+   return exp(log(tbase(a)) * b);
+}
+
+template<>
+inline FDual<double> pow(const double &a, const FDual<double> &b)
+{
+   return exp(std::log(a) * b);
+}
+
+template<typename tbase>
+inline FDual<tbase> sin(const FDual<tbase> &f)
+{
+   return FDual<tbase>(sin(f.real()), f.dual() * cos(f.real()));
+}
+
+template<>
+inline FDual<double> sin(const FDual<double> &f)
+{
+   return FDual<double>(std::sin(f.real()), f.dual() * std::cos(f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> sinh(const FDual<tbase> &f)
+{
+   return FDual<tbase>(sinh(f.real()), f.dual() * cosh(f.real()));
+}
+
+template<>
+inline FDual<double> sinh(const FDual<double> &f)
+{
+   return FDual<double>(std::sinh(f.real()), f.dual() * std::cosh(f.real()));
+}
+
+template<typename tbase>
+inline FDual<tbase> sqrt(const FDual<tbase> &f)
+{
+   tbase a = sqrt(f.real());
+   return FDual<tbase>(a, f.dual() / (tbase(2) * a));
+}
+
+template<>
+inline FDual<double> sqrt(const FDual<double> &f)
+{
+   double a = std::sqrt(f.real());
+   return FDual<double>(a, f.dual() / (double(2) * a));
+}
+
+template<typename tbase>
+inline FDual<tbase> tan(const FDual<tbase> &f)
+{
+   tbase a = tan(f.real());
+   return FDual<tbase>(a, f.dual() * (tbase(1) + a * a));
+}
+
+template<>
+inline FDual<double> tan(const FDual<double> &f)
+{
+   double a = std::tan(f.real());
+   return FDual<double>(a, f.dual() * (double(1) + a * a));
+}
+
+template<typename tbase>
+inline FDual<tbase> tanh(const FDual<tbase> &f)
+{
+   tbase a = tanh(f.real());
+   return FDual<tbase>(a, f.dual() * (tbase(1) - a * a));
+}
+
+template<>
+inline FDual<double> tanh(const FDual<double> &f)
+{
+   double a = std::tanh(f.real());
+   return FDual<double>(a, f.dual() * (double(1) - a * a));
+}
+
+} // namespace ad
+
+} // namespace mfem
+
+#endif

linalg/taddensemat.hpp

@@ -0,0 +1,516 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef TADDENSEMATRIX_H
+#define TADDENSEMATRIX_H
+
+#include "../config/config.hpp"
+#include "../general/globals.hpp"
+
+#include "tadvector.hpp"
+#include "densemat.hpp"
+
+namespace mfem
+{
+/// Templated dense matrix data type.
+/** The main goal of the TADDenseMatrix class is to serve as a data
+  container for representing dense matrices in classes, methods, and
+  functions utilized with automatic differentiation (AD). The
+  functionality/interface is copied from the standard MFEM dense
+  matrix mfem::DenseMatrix. The basic idea is to utilize the templated
+  vector class in combination with AD during the development phase.
+  The AD parts can be replaced with optimized code once the initial
+  development of the application is complete.  The common interface
+  between TADDenseMatrix and DenseMatrix will ease the transition
+  from AD to hand-optimized code as it does not require a change
+  in the interface or the code structure. TADDenseMatrix is intended
+  to be utilized for dense serial matrices. The objects can be combined
+  with TADVector or standard Vector.*/
+template<typename dtype>
+class TADDenseMatrix
+{
+private:
+   int height; ///< Dimension of the output / number of rows in the matrix.
+   int width;  ///< Dimension of the input / number of columns in the matrix.
+   dtype *data;
+   int capacity; // zero or negative capacity means we do not own the data.
+
+public:
+   /// Get the height (size of output) of the Operator. Synonym with NumRows().
+   inline int Height() const { return height; }
+   /** @brief Get the number of rows (size of output) of the Operator. Synonym
+       with Height(). */
+   inline int NumRows() const { return height; }
+
+   /// Get the width (size of input) of the Operator. Synonym with NumCols().
+   inline int Width() const { return width; }
+   /** @brief Get the number of columns (size of input) of the Operator. Synonym
+       with Width(). */
+   inline int NumCols() const { return width; }
+
+   /** Default constructor for TADDenseMatrix.
+      Sets data = NULL and height = width = 0. */
+   TADDenseMatrix()
+   {
+      data = nullptr;
+      capacity = 0;
+      height = 0;
+      width = 0;
+   }
+
+   /// Copy constructor
+   template<typename idtype>
+   TADDenseMatrix(const TADDenseMatrix<idtype> &m)
+   {
+      height = m.GetHeight();
+      width = m.GetWidth();
+      const int hw = height * width;
+      if (hw > 0)
+      {
+         idtype *mdata = m.Data();
+         MFEM_ASSERT(mdata, "invalid source matrix");
+         data = new dtype[hw];
+         capacity = hw;
+         for (int i = 0; i < hw; i++)
+         {
+            data[i] = mdata[i];
+         }
+      }
+      else
+      {
+         data = nullptr;
+         capacity = 0;
+         width = 0;
+         height = 0;
+      }
+   }
+   /// Copy constructor using standard DenseMatrix
+   TADDenseMatrix(const DenseMatrix &m)
+   {
+      height = m.Height();
+      width = m.Width();
+      const int hw = height * width;
+      if (hw > 0)
+      {
+         double *mdata = m.Data();
+         MFEM_ASSERT(mdata, "invalid source matrix");
+         data = new dtype[hw];
+         capacity = hw;
+         for (int i = 0; i < hw; i++)
+         {
+            data[i] = mdata[i];
+         }
+      }
+      else
+      {
+         data = nullptr;
+         capacity = 0;
+         width = 0;
+         height = 0;
+      }
+   }
+
+   /// Creates square matrix of size s.
+   explicit TADDenseMatrix(int s)
+   {
+      MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
+      height = s;
+      width = s;
+      capacity = s * s;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity](); // init with zeroes
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   /// Creates rectangular matrix of size m x n.
+   TADDenseMatrix(int m, int n)
+   {
+      MFEM_ASSERT(m >= 0 && n >= 0,
+                  "invalid DenseMatrix size: " << m << " x " << n);
+      height = m;
+      width = n;
+      capacity = m * n;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity](); // init with zeroes
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   TADDenseMatrix(const TADDenseMatrix<dtype> &mat, char ch)
+   {
+      height = mat.Width();
+      width = mat.Height();
+      capacity = height * width;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity];
+
+         for (int i = 0; i < height; i++)
+         {
+            for (int j = 0; j < width; j++)
+            {
+               (*this)(i, j) = mat(j, i);
+            }
+         }
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   /// Change the size of the DenseMatrix to s x s.
+   void SetSize(int s) { SetSize(s, s); }
+
+   /// Change the size of the DenseMatrix to h x w.
+   void SetSize(int h, int w)
+   {
+      MFEM_ASSERT(h >= 0 && w >= 0,
+                  "invalid DenseMatrix size: " << h << " x " << w);
+      if (Height() == h && Width() == w)
+      {
+         return;
+      }
+      height = h;
+      width = w;
+      const int hw = h * w;
+      if (hw > std::abs(capacity))
+      {
+         if (capacity > 0)
+         {
+            delete[] data;
+         }
+         capacity = hw;
+         data = new dtype[hw](); // init with zeroes
+      }
+   }
+
+   /// Returns the matrix data array.
+   inline dtype *Data() const { return data; }
+   /// Returns the matrix data array.
+   inline dtype *GetData() const { return data; }
+
+   inline bool OwnsData() const { return (capacity > 0); }
+
+   /// Returns reference to a_{ij}.
+   dtype &operator()(int i, int j)
+   {
+      MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+      return data[i + j * height];
+   }
+
+   const dtype &operator()(int i, int j) const
+   {
+      MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+      return data[i + j * height];
+   }
+
+   dtype &Elem(int i, int j) { return (*this)(i, j); }
+
+   const dtype &Elem(int i, int j) const { return (*this)(i, j); }
+
+   void Mult(const dtype *x, dtype *y) const
+   {
+      if (width == 0)
+      {
+         for (int row = 0; row < height; row++)
+         {
+            y[row] = 0.0;
+         }
+         return;
+      }
+      dtype *d_col = data;
+      dtype x_col = x[0];
+      for (int row = 0; row < height; row++)
+      {
+         y[row] = x_col * d_col[row];
+      }
+      d_col += height;
+      for (int col = 1; col < width; col++)
+      {
+         x_col = x[col];
+         for (int row = 0; row < height; row++)
+         {
+            y[row] += x_col * d_col[row];
+         }
+         d_col += height;
+      }
+   }
+
+   void Mult(const TADVector<dtype> &x, TADVector<dtype> &y) const
+   {
+      MFEM_ASSERT(height == y.Size() && width == x.Size(),
+                  "incompatible dimensions");
+
+      Mult((const dtype *) x, (dtype *) y);
+   }
+
+   dtype operator*(const TADDenseMatrix<dtype> &m) const
+   {
+      MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
+                  "incompatible dimensions");
+
+      const int hw = height * width;
+      dtype a = 0.0;
+      for (int i = 0; i < hw; i++)
+      {
+         a += data[i] * m.data[i];
+      }
+
+      return a;
+   }
+
+   void MultTranspose(const dtype *x, dtype *y) const
+   {
+      dtype *d_col = data;
+      for (int col = 0; col < width; col++)
+      {
+         double y_col = 0.0;
+         for (int row = 0; row < height; row++)
+         {
+            y_col += x[row] * d_col[row];
+         }
+         y[col] = y_col;
+         d_col += height;
+      }
+   }
+
+   void MultTranspose(const TADVector<dtype> &x, TADVector<dtype> &y) const
+   {
+      MFEM_ASSERT(height == x.Size() && width == y.Size(),
+                  "incompatible dimensions");
+
+      MultTranspose((const dtype *) x, (dtype *) y);
+   }
+
+   void Randomize(int seed)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double) (RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int) time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned) seed);
+
+      for (int i = 0; i < capacity; i++)
+      {
+         data[i] = (dtype)(std::abs(rand() / max));
+      }
+   }
+
+   void RandomizeDiag(int seed)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double) (RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int) time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned) seed);
+
+      for (int i = 0; i < std::min(height, width); i++)
+      {
+         Elem(i, i) = (dtype)(std::abs(rand() / max));
+      }
+   }
+
+   /// Creates n x n diagonal matrix with diagonal elements c
+   void Diag(dtype c, int n)
+   {
+      SetSize(n);
+
+      const int N = n * n;
+      for (int i = 0; i < N; i++)
+      {
+         data[i] = (dtype) 0.0;
+      }
+      for (int i = 0; i < n; i++)
+      {
+         data[i * (n + 1)] = c;
+      }
+   }
+   /// Creates n x n diagonal matrix with diagonal given by diag
+   template<typename itype>
+   void Diag(itype *diag, int n)
+   {
+      SetSize(n);
+
+      int i, N = n * n;
+      for (i = 0; i < N; i++)
+      {
+         data[i] = 0.0;
+      }
+      for (i = 0; i < n; i++)
+      {
+         data[i * (n + 1)] = (dtype) diag[i];
+      }
+   }
+
+   /// (*this) = (*this)^t
+   void Transpose()
+   {
+      int i, j;
+      dtype t;
+
+      if (Width() == Height())
+      {
+         for (i = 0; i < Height(); i++)
+            for (j = i + 1; j < Width(); j++)
+            {
+               t = (*this)(i, j);
+               (*this)(i, j) = (*this)(j, i);
+               (*this)(j, i) = t;
+            }
+      }
+      else
+      {
+         TADDenseMatrix<dtype> T(*this, 't');
+         (*this) = T;
+      }
+   }
+   /// (*this) = A^t
+   template<typename itype>
+   void Transpose(const TADDenseMatrix<itype> &A)
+   {
+      SetSize(A.Width(), A.Height());
+
+      for (int i = 0; i < Height(); i++)
+         for (int j = 0; j < Width(); j++)
+         {
+            (*this)(i, j) = (dtype) A(j, i);
+         }
+   }
+
+   /// (*this) = 1/2 ((*this) + (*this)^t)
+   void Symmetrize()
+   {
+#ifdef MFEM_DEBUG
+      if (Width() != Height())
+      {
+         mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
+      }
+#endif
+
+      for (int i = 0; i < Height(); i++)
+         for (int j = 0; j < i; j++)
+         {
+            dtype a = 0.5 * ((*this)(i, j) + (*this)(j, i));
+            (*this)(j, i) = (*this)(i, j) = a;
+         }
+   }
+
+   void Lump()
+   {
+      for (int i = 0; i < Height(); i++)
+      {
+         dtype L = 0.0;
+         for (int j = 0; j < Width(); j++)
+         {
+            L += (*this)(i, j);
+            (*this)(i, j) = (dtype) 0.0;
+         }
+         (*this)(i, i) = L;
+      }
+   }
+};
+
+template<typename dtype>
+void CalcAdjugate(const TADDenseMatrix<dtype> &a, TADDenseMatrix<dtype> &adja)
+{
+#ifdef MFEM_DEBUG
+   if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
+   {
+      mfem_error("CalcAdjugate(...)");
+   }
+   if (a.Width() != adja.Height() || a.Height() != adja.Width())
+   {
+      mfem_error("CalcAdjugate(...)");
+   }
+#endif
+
+   if (a.Width() < a.Height())
+   {
+      const dtype *d = a.Data();
+      dtype *ad = adja.Data();
+      if (a.Width() == 1)
+      {
+         // N x 1, N = 2,3
+         ad[0] = d[0];
+         ad[1] = d[1];
+         if (a.Height() == 3)
+         {
+            ad[2] = d[2];
+         }
+      }
+      else
+      {
+         // 3 x 2
+         double e, g, f;
+         e = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
+         g = d[3] * d[3] + d[4] * d[4] + d[5] * d[5];
+         f = d[0] * d[3] + d[1] * d[4] + d[2] * d[5];
+
+         ad[0] = d[0] * g - d[3] * f;
+         ad[1] = d[3] * e - d[0] * f;
+         ad[2] = d[1] * g - d[4] * f;
+         ad[3] = d[4] * e - d[1] * f;
+         ad[4] = d[2] * g - d[5] * f;
+         ad[5] = d[5] * e - d[2] * f;
+      }
+      return;
+   }
+
+   if (a.Width() == 1)
+   {
+      adja(0, 0) = (dtype) 1.0;
+   }
+   else if (a.Width() == 2)
+   {
+      adja(0, 0) = a(1, 1);
+      adja(0, 1) = -a(0, 1);
+      adja(1, 0) = -a(1, 0);
+      adja(1, 1) = a(0, 0);
+   }
+   else
+   {
+      adja(0, 0) = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1);
+      adja(0, 1) = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2);
+      adja(0, 2) = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1);
+
+      adja(1, 0) = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2);
+      adja(1, 1) = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0);
+      adja(1, 2) = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2);
+
+      adja(2, 0) = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0);
+      adja(2, 1) = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1);
+      adja(2, 2) = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
+   }
+}
+
+} // namespace mfem
+
+#endif

linalg/tadvector.hpp

@@ -0,0 +1,700 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_TADVECTOR
+#define MFEM_TADVECTOR
+
+#include "../general/mem_manager.hpp"
+#include "vector.hpp"
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+#if defined(_MSC_VER) && (_MSC_VER < 1800)
+#include <float.h>
+#define isfinite _finite
+#endif
+
+namespace mfem
+{
+/// Templated vector data type.
+/** The main goal of the TADVector class is to serve as a data
+  container for representing vectors in classes, methods, and
+  functions utilized with automatic differentiation (AD). The
+  functionality/interface is copied from the standard MFEM dense
+  vector mfem::Vector.  The basic idea is to utilize the templated
+  vector class in combination with AD during the development phase.
+  The AD parts can be replaced with optimized code once the initial
+  development of the application is complete.  The common interface
+  between TADVector and Vector will ease the transition from AD to
+  hand-optimized code as it does not require a change in the
+  interface or the code structure. TADVector is intended to be
+  utilized for dense serial vectors. */
+template<typename dtype>
+class TADVector
+{
+protected:
+   dtype *data;
+   int size;
+   int capacity;
+
+public:
+   /// Default constructor for Vector. Sets size = 0 and data = NULL.
+   TADVector()
+   {
+      data = nullptr;
+      size = 0;
+      capacity = 0;
+   }
+
+   /// Copy constructor. Allocates a new data array and copies the data.
+   TADVector(const TADVector<dtype> &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         data = new dtype[s];
+         capacity = s;
+         for (int i = 0; i < s; i++)
+         {
+            data[i] = v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         capacity = 0;
+         data = nullptr;
+      }
+   }
+
+   TADVector(const Vector &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         capacity = s;
+         data = new dtype[s];
+         for (int i = 0; i < s; i++)
+         {
+            data[i] = v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         capacity = 0;
+         data = nullptr;
+      }
+   }
+
+   /// @brief Creates vector of size s.
+   /// @warning Entries are not initialized to zero!
+   explicit TADVector(int s)
+   {
+      if (s > 0)
+      {
+         size = s;
+         capacity = s;
+         data = new dtype[size];
+      }
+      else
+      {
+         size = 0;
+         capacity = 0;
+         data = nullptr;
+      }
+   }
+
+   /// Creates a vector referencing an array of doubles, owned by someone else.
+   /** The pointer @a _data can be NULL. The data array can be replaced later
+       with SetData(). */
+   TADVector(dtype *_data, int _size)
+   {
+      if (capacity > 0)
+      {
+         delete[] data;
+         capacity = 0;
+      }
+      size = _size;
+      data = _data;
+   }
+
+   /// Reads a vector from multiple files
+   void Load(std::istream **in, int np, int *dim)
+   {
+      int i, j, s;
+
+      s = 0;
+      for (i = 0; i < np; i++)
+      {
+         s += dim[i];
+      }
+      SetSize(s);
+
+      int p = 0;
+      double tmpd;
+      for (i = 0; i < np; i++)
+      {
+         for (j = 0; j < dim[i]; j++)
+         {
+            *in[i] >> tmpd;
+            data[p++] = dtype(tmpd);
+         }
+      }
+   }
+
+   /// Load a vector from an input stream.
+   void Load(std::istream &in, int Size)
+   {
+      SetSize(Size);
+      double tmpd;
+      for (int i = 0; i < size; i++)
+      {
+         in >> tmpd;
+         data[i] = dtype(tmpd);
+      }
+   }
+
+   /// Load a vector from an input stream, reading the size from the stream.
+   void Load(std::istream &in)
+   {
+      int s;
+      in >> s;
+      Load(in, s);
+   }
+
+   /// @brief Resize the vector to size @a s.
+   /** If the new size is less than or equal to Capacity() then the internal
+       data array remains the same. Otherwise, the old array is deleted, if
+       owned, and a new array of size @a s is allocated without copying the
+       previous content of the Vector.
+       @warning In the second case above (new size greater than current one),
+       the vector will allocate new data array, even if it did not own the
+       original data! Also, new entries are not initialized! */
+   void SetSize(int s)
+   {
+      if (s == size)
+      {
+         return;
+      }
+
+      if (s <= capacity)
+      {
+         size = s;
+         return;
+      }
+
+      delete[] data;
+      data = new dtype[s];
+      size = s;
+      capacity = s;
+   }
+
+   /// Set the Vector data and size.
+   /** The Vector does not assume ownership of the new data. The new size is
+       @warning This method should be called only when OwnsData() is false.
+       @sa NewDataAndSize(). */
+   void SetDataAndSize(dtype *d, int s)
+   {
+      if (OwnsData())
+      {
+         delete[] data;
+         capacity = 0;
+      }
+      data = d;
+      size = s;
+   }
+
+   /// Set the Vector data and size, deleting the old data, if owned.
+   /** The Vector does not assume ownership of the new data. The new size is
+       also used as the new Capacity().
+       @sa SetDataAndSize(). */
+   void NewDataAndSize(dtype *d, int s) { SetDataAndSize(d, s); }
+
+   /// Reset the Vector to be a reference to a sub-vector of @a base.
+   inline void MakeRef(TADVector<dtype> &base, int offset, int size_)
+   {
+      NewDataAndSize(base.GetData() + offset, size_);
+   }
+
+   /** @brief Reset the Vector to be a reference to a sub-vector of @a base
+       without changing its current size. */
+   inline void MakeRef(TADVector<dtype> &base, int offset)
+   {
+      int tsiz = size;
+      NewDataAndSize(base.GetData() + offset, tsiz);
+   }
+
+   /// Destroy a vector
+   void Destroy()
+   {
+      size = 0;
+      capacity = 0;
+      delete[] data;
+   }
+
+   /// Returns the size of the vector.
+   inline int Size() const { return size; }
+
+   /// Return the size of the currently allocated data array.
+   /** It is always true that Capacity() >= Size(). */
+   inline int Capacity() const { return capacity; }
+
+   /// Return a pointer to the beginning of the Vector data.
+   /** @warning This method should be used with caution as it gives write access
+       to the data of const-qualified Vector%s. */
+   inline dtype *GetData() const
+   {
+      return const_cast<dtype *>((const dtype *) data);
+   }
+
+   /// Conversion to `double *`.
+   /** @note This conversion function makes it possible to use [] for indexing
+       in addition to the overloaded operator()(int). */
+   inline operator dtype *() { return data; }
+
+   /// Conversion to `const double *`.
+   /** @note This conversion function makes it possible to use [] for indexing
+       in addition to the overloaded operator()(int). */
+   inline operator const dtype *() const { return data; }
+
+   /// Read the Vector data (host pointer) ownership flag.
+   inline bool OwnsData() const { return (capacity > 0); }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline void StealData(dtype **p)
+   {
+      *p = data;
+      delete[] data;
+      size = 0;
+      capacity = 0;
+   }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline dtype *StealData()
+   {
+      dtype *p;
+      StealData(&p);
+      return p;
+   }
+
+   /// Access Vector entries. Index i = 0 .. size-1.
+   dtype &Elem(int i) { return operator()(i); }
+   /// Read only access to Vector entries. Index i = 0 .. size-1.
+   const dtype &Elem(int i) const { return operator()(i); }
+
+   /// Access Vector entries using () for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline dtype &operator()(int i)
+   {
+      MFEM_ASSERT(data && i >= 0 && i < size,
+                  "index [" << i << "] is out of range [0," << size << ")");
+
+      return data[i];
+   }
+
+   /// Read only access to Vector entries using () for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline const dtype &operator()(int i) const
+   {
+      MFEM_ASSERT(data && i >= 0 && i < size,
+                  "index [" << i << "] is out of range [0," << size << ")");
+
+      return data[i];
+   }
+
+   /// Dot product with a `dtype *` array.
+   dtype operator*(const dtype *v) const
+   {
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   /// Return the inner-product.
+   dtype operator*(const TADVector<dtype> &v) const
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   dtype operator*(const Vector &v) const
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   /// Copy Size() entries from @a v.
+   TADVector<dtype> &operator=(const dtype *v)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = v[i];
+      }
+      return *this;
+   }
+
+   /// Copy assignment.
+   /** @note Defining this method overwrites the implicitly defined copy
+       assignemnt operator. */
+   TADVector<dtype> &operator=(const TADVector<dtype> &v)
+   {
+      SetSize(v.Size());
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = v[i];
+      }
+      return *this;
+   }
+
+   TADVector<dtype> &operator=(const Vector &v)
+   {
+      SetSize(v.Size());
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = v[i];
+      }
+      return *this;
+   }
+
+   /// Redefine '=' for vector = constant.
+   template<typename ivtype>
+   TADVector &operator=(ivtype value)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = (dtype) value;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator*=(ivtype c)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] * c;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator/=(ivtype c)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] / c;
+      }
+      return *this;
+   }
+
+   TADVector &operator-=(const TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] - v[i];
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator-=(ivtype v)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] - v;
+      }
+      return *this;
+   }
+
+   TADVector &operator+=(const TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] + v[i];
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator+=(ivtype v)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] + v;
+      }
+      return *this;
+   }
+
+   /// (*this) += a * Va
+   template<typename ivtype, typename vtype>
+   TADVector &Add(const ivtype a, const vtype &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] + a * v[i];
+      }
+      return *this;
+   }
+
+   /// (*this) = a * x
+   template<typename ivtype, typename vtype>
+   TADVector &Set(const ivtype a, const vtype &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = a * v[i];
+      }
+      return *this;
+   }
+
+   template<typename vtype>
+   void SetVector(const vtype &v, int offset)
+   {
+      MFEM_ASSERT(v.Size() + offset <= size, "invalid sub-vector");
+      for (int i = 0; i < size; i++)
+      {
+         data[i + offset] = v[i];
+      }
+   }
+
+   /// (*this) = -(*this)
+   void Neg()
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = -data[i];
+      }
+   }
+
+   /// Swap the contents of two Vectors
+   inline void Swap(TADVector<dtype> &other)
+   {
+      Swap(data, other.data);
+      Swap(size, other.size);
+      Swap(capacity, other.capacity);
+   }
+
+   /// Set v = v1 + v2.
+   template<typename vtype1, typename vtype2>
+   friend void add(const vtype1 &v1, const vtype2 &v2, TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(v1.Size() == v.Size(), "incompatible Vectors!");
+      MFEM_ASSERT(v2.Size() == v.Size(), "incompatible Vectors!");
+      for (int i = 0; i < v.Size(); i++)
+      {
+         v[i] = v1[i] + v2[i];
+      }
+   }
+
+   /// Set v = v1 + alpha * v2.
+   template<typename vtype1, typename vtype2>
+   friend void add(const vtype1 &v1,
+                   dtype alpha,
+                   const vtype2 &v2,
+                   TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(v1.Size() == v.Size(), "incompatible Vectors!");
+      MFEM_ASSERT(v2.Size() == v.Size(), "incompatible Vectors!");
+      for (int i = 0; i < v.Size(); i++)
+      {
+         v[i] = v1[i] + alpha * v2[i];
+      }
+   }
+
+   template<typename vtype1, typename vtype2>
+   friend void add(const dtype a,
+                   const vtype1 &x,
+                   const dtype b,
+                   const vtype2 &y,
+                   TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = a * x[i] + b * y[i];
+      }
+   }
+
+   template<typename vtype1, typename vtype2>
+   friend void add(const dtype a,
+                   const vtype1 &x,
+                   const vtype2 &y,
+                   TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = a * x[i] + y[i];
+      }
+   }
+
+   template<typename vtype1, typename vtype2>
+   friend void subtract(const vtype1 &x, const vtype2 &y, TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = x[i] - y[i];
+      }
+   }
+
+   template<typename ivtype, typename vtype1, typename vtype2>
+   friend void subtract(const ivtype a,
+                        const vtype1 &x,
+                        const vtype2 &y,
+                        TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = a * (x[i] - y[i]);
+      }
+   }
+
+   /// Destroys vector.
+   ~TADVector() { delete[] data; }
+
+   /// Prints vector to stream out.
+   void Print(std::ostream &out = mfem::out, int width = 8) const
+   {
+      if (!size)
+      {
+         return;
+      }
+      for (int i = 0; 1;)
+      {
+         out << data[i];
+         i++;
+         if (i == size)
+         {
+            break;
+         }
+         if (i % width == 0)
+         {
+            out << '\n';
+         }
+         else
+         {
+            out << ' ';
+         }
+      }
+      out << '\n';
+   }
+
+   /// Set random values in the vector.
+   void Randomize(int seed = 0)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double) (RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int) time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned) seed);
+
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = std::abs(rand() / max);
+      }
+   }
+   /// Returns the l2 norm of the vector.
+   dtype Norml2() const
+   {
+      // Scale entries of Vector on the fly, using algorithms from
+      // std::hypot() and LAPACK's drm2. This scaling ensures that the
+      // argument of each call to std::pow is <= 1 to avoid overflow.
+      if (0 == size)
+      {
+         return 0.0;
+      } // end if 0 == size
+
+      if (1 == size)
+      {
+         return abs(data[0]);
+      } // end if 1 == size
+
+      dtype scale = 0.0;
+      dtype sum = 0.0;
+
+      for (int i = 0; i < size; i++)
+      {
+         if (data[i] != 0.0)
+         {
+            const dtype absdata = abs(data[i]);
+            if (scale <= absdata)
+            {
+               const dtype sqr_arg = scale / absdata;
+               sum = 1.0 + sum * (sqr_arg * sqr_arg);
+               scale = absdata;
+               continue;
+            } // end if scale <= absdata
+            const dtype sqr_arg = absdata / scale;
+            sum += (sqr_arg * sqr_arg); // else scale > absdata
+         }                              // end if data[i] != 0
+      }
+      return scale * sqrt(sum);
+   }
+
+   /// Returns the l_infinity norm of the vector.
+   dtype Normlinf() const
+   {
+      dtype max = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         max = max(abs(data[i]), max);
+      }
+      return max;
+   }
+   /// Returns the l_1 norm of the vector.
+   dtype Norml1() const
+   {
+      dtype sum = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         sum += abs(data[i]);
+      }
+      return sum;
+   }
+};
+
+} // namespace mfem
+
+#endif

makefile

@@ -267,7 +267,8 @@ endif
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
- OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX
+ OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX ADEPT FADBADPP
+
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
@@ -333,7 +334,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_OCCA\
  MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2\
  MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_USE_MUMPS MFEM_SOURCE_DIR\
- MFEM_INSTALL_DIR
+ MFEM_INSTALL_DIR MFEM_USE_ADEPT MFEM_USE_FADBADPP MFEM_USE_ADFORWARD
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -657,6 +658,9 @@ status info:
 	$(info MFEM_USE_UMPIRE        = $(MFEM_USE_UMPIRE))
 	$(info MFEM_USE_SIMD          = $(MFEM_USE_SIMD))
 	$(info MFEM_USE_ADIOS2        = $(MFEM_USE_ADIOS2))
+	$(info MFEM_USE_ADEPT         = $(MFEM_USE_ADEPT))
+	$(info MFEM_USE_FADBADPP      = $(MFEM_USE_FADBADPP))
+	$(info MFEM_USE_ADFORWARD     = $(MFEM_USE_ADFORWARD))
 	$(info MFEM_USE_MKL_CPARDISO  = $(MFEM_USE_MKL_CPARDISO))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_HOST_CXX          = $(value MFEM_HOST_CXX))

tests/unit/CMakeLists.txt

@@ -31,6 +31,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_matrix_sparse.cpp
   linalg/test_matrix_square.cpp
   linalg/test_ode.cpp
+  linalg/test_fdual.cpp
   linalg/test_ode2.cpp
   linalg/test_operator.cpp
   linalg/test_cg_indefinite.cpp

tests/unit/linalg/test_fdual.cpp

@@ -0,0 +1,203 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+using namespace mfem;
+
+template<typename tbase>
+tbase exprp02(tbase x, tbase y)
+{
+   return sin(x) * cos(y) + tan(x * y);
+}
+
+template<typename tbase>
+tbase exprp02x(tbase x, tbase y)
+{
+   return cos(x) * cos(y) + y * (1.0 + pow(tan(x * y), 2.0));
+}
+
+template<typename tbase>
+tbase exprp02y(tbase x, tbase y)
+{
+   return -sin(x) * sin(y) + x * (1.0 + pow(tan(x * y), 2.0));
+}
+
+TEST_CASE("Simple AD tests", "[Simple_AD_tests]")
+{
+   SECTION("sin")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::sin(xx);
+      d = std::cos(x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("cos")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::cos(xx);
+      d = -std::sin(x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("tan")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::tan(xx);
+      d = 1.0 + std::tan(x) * std::tan(x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("exp")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::exp(xx);
+      d = exp(x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("log")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::log(xx);
+      d = 1.0 / x;
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("pow")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::pow(xx, 1.5);
+      d = 1.5 * std::pow(x, 0.5);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("atan")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::atan(xx);
+      d = 1.0 / (1.0 + x * x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("asin")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::asin(xx);
+      d = 1.0 / std::sqrt(1.0 - x * x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("acos")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x, 1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::acos(xx);
+      d = -1.0 / std::sqrt(1.0 - x * x);
+      REQUIRE(std::abs(d - lrez.dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("general")
+   {
+      double x = 1.0;
+      double y = 1.5;
+
+      double pr = exprp02(x, y);
+      double dx = exprp02x(x, y);
+      double dy = exprp02y(x, y);
+
+      {
+         ad::FDual<double> xx(x, 1.0);
+         ad::FDual<double> yy(y, 0.0);
+         ad::FDual<double> rr = exprp02(xx, yy);
+         REQUIRE(std::abs(rr.real() - pr)
+                 < std::numeric_limits<double>::epsilon());
+         REQUIRE(std::abs(rr.dual() - dx)
+                 < std::numeric_limits<double>::epsilon());
+      }
+
+      {
+         ad::FDual<double> xx(x, 0.0);
+         ad::FDual<double> yy(y, 1.0);
+         ad::FDual<double> rr = exprp02(xx, yy);
+         REQUIRE(std::abs(rr.real() - pr)
+                 < std::numeric_limits<double>::epsilon());
+         REQUIRE(std::abs(rr.dual() - dy)
+                 < std::numeric_limits<double>::epsilon());
+      }
+   }
+
+   SECTION("second_derivative")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<ad::FDual<double>> xxx(ad::FDual<double>(x,
+                                                         1.0),
+                                       ad::FDual<double>(1.0,
+                                                         0.0));
+      ad::FDual<ad::FDual<double>> drez = ad::exp(xxx);
+      d = exp(x);
+      REQUIRE(std::abs(d - drez.dual().dual())
+              < std::numeric_limits<double>::epsilon());
+
+      drez = ad::log(xxx);
+      d = -1.0 / (x * x);
+      REQUIRE(std::abs(d - drez.dual().dual())
+              < std::numeric_limits<double>::epsilon());
+
+      drez = ad::sin(xxx);
+      d = -sin(x);
+      REQUIRE(std::abs(d - drez.dual().dual())
+              < std::numeric_limits<double>::epsilon());
+
+      drez = ad::cos(xxx);
+      d = -cos(x);
+      REQUIRE(std::abs(d - drez.dual().dual())
+              < std::numeric_limits<double>::epsilon());
+   }
+}

tests/unit/linalg/test_vector.cpp

@@ -123,3 +123,129 @@ TEST_CASE("Vector Tests", "[Vector]")
       REQUIRE(diff.Norml2() < tol);
    }
 }
+
+TEST_CASE("TADVector Tests", "[TADVector]")
+{
+   double tol = 1e-12;
+
+   TADVector<double> a(3), b(3);
+   Vector bm(3);
+   bm(0) = 2.0;
+   bm(1) = 1.0;
+   bm(2) = 4.0;
+
+   a(0) = 1.0;
+   a(1) = 3.0;
+   a(2) = 5.0;
+
+   b(0) = 2.0;
+   b(1) = 1.0;
+   b(2) = 4.0;
+
+   double bp[3];
+   bp[0] = b(0);
+   bp[1] = b(1);
+   bp[2] = b(2);
+
+   TADVector<double> apb(3), amb(3);
+   Vector amm(3);
+   apb(0) = 3.0;
+   apb(1) = 4.0;
+   apb(2) = 9.0;
+
+   amb(0) = -1.0;
+   amb(1) = 2.0;
+   amb(2) = 1.0;
+
+   amm(0) = -1.0;
+   amm(1) = 2.0;
+   amm(2) = 1.0;
+
+   TADVector<double> tmp(3);
+   TADVector<double> diff(3);
+
+   SECTION("Dot product")
+   {
+      REQUIRE(a * b - 25.0 < tol);
+      REQUIRE(a * bp - 25.0 < tol);
+      REQUIRE(a * bm - 25.0 < tol);
+   }
+
+   SECTION("Multiply and divide")
+   {
+      a *= 3.0;
+      b /= -4.0;
+
+      REQUIRE(a * b + 3.0 * 25.0 / 4.0 < tol);
+      REQUIRE(a * bp - 3.0 * 25.0 < tol);
+      REQUIRE(a * bm - 3.0 * 25.0 < tol);
+   }
+
+   SECTION("Minus scalar")
+   {
+      a -= 3.0;
+      REQUIRE(a.Norml2() - sqrt(8.0) < tol);
+   }
+
+   SECTION("Minus vector")
+   {
+      a -= b;
+      subtract(a, amb, diff);
+      REQUIRE(diff.Norml2() < tol);
+      subtract(a, amm, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+   SECTION("Subtract vector")
+   {
+      subtract(0.5, a, b, tmp);
+      tmp *= 2.0;
+      subtract(tmp, amb, diff);
+      REQUIRE(diff.Norml2() < tol);
+      subtract(tmp, amm, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+   SECTION("Plus scalar")
+   {
+      a += 2.0;
+      REQUIRE(a.Norml2() - sqrt(83.0) < tol);
+   }
+
+   SECTION("Plus vector")
+   {
+      a += b;
+      add(1.0, a, -1.0, apb, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+   SECTION("Add vector 1")
+   {
+      a.Add(1.0, b);
+      add(1.0, a, -1.0, apb, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+   SECTION("Add vector 2")
+   {
+      add(a, b, tmp);
+      apb.Neg();
+      add(tmp, apb, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+   SECTION("Add vector 3")
+   {
+      add(a, 1.0, b, tmp);
+      apb.Neg();
+      add(tmp, apb, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+   SECTION("Add vector 4")
+   {
+      add(1.0, a, b, tmp);
+      subtract(tmp, apb, diff);
+      REQUIRE(diff.Norml2() < tol);
+   }
+}