45template <
typename ScalarT>
46ScalarT
func(
const ScalarT&
a,
const ScalarT& b,
const ScalarT&
c) {
47 ScalarT r =
c*std::log(b+1.)/std::sin(
a);
52void func_deriv(
double a,
double b,
double c,
double& drda,
double& drdb)
54 drda = -(
c*std::log(b+1.)/std::pow(std::sin(
a),2.))*std::cos(
a);
55 drdb =
c / ((b+1.)*std::sin(
a));
59void func_deriv2(
double a,
double b,
double c,
double& d2rda2,
double& d2rdb2,
62 d2rda2 =
c*std::log(b+1.)/std::sin(
a) + 2.*(
c*std::log(b+1.)/std::pow(std::sin(
a),3.))*std::pow(std::cos(
a),2.);
63 d2rdb2 = -
c / (std::pow(b+1.,2.)*std::sin(
a));
64 d2rdadb = -
c / ((b+1.)*std::pow(std::sin(
a),2.))*std::cos(
a);
67int main(
int argc,
char **argv)
69 double pi = std::atan(1.0)*4.0;
97 double d2rda2, d2rdb2, d2rdadb;
101 rfad =
func(afad, bfad, cfad);
104 double r_ad = rfad.val().val();
105 double drda_ad = rfad.dx(0).val();
106 double drdb_ad = rfad.dx(1).val();
107 double d2rda2_ad = rfad.dx(0).dx(0);
108 double d2rdadb_ad = rfad.dx(0).dx(1);
109 double d2rdbda_ad = rfad.dx(1).dx(0);
110 double d2rdb2_ad = rfad.dx(1).dx(1);
115 std::cout.setf(std::ios::scientific);
116 std::cout.precision(
p);
117 std::cout <<
" r = " << std::setw(w) << r <<
" (original) == "
118 << std::setw(w) << r_ad <<
" (AD) Error = " << std::setw(w)
119 << r - r_ad << std::endl
120 <<
" dr/da = " << std::setw(w) << drda <<
" (analytic) == "
121 << std::setw(w) << drda_ad <<
" (AD) Error = " << std::setw(w)
122 << drda - drda_ad << std::endl
123 <<
" dr/db = " << std::setw(w) << drdb <<
" (analytic) == "
124 << std::setw(w) << drdb_ad <<
" (AD) Error = " << std::setw(w)
125 << drdb - drdb_ad << std::endl
126 <<
"d^2r/da^2 = " << std::setw(w) << d2rda2 <<
" (analytic) == "
127 << std::setw(w) << d2rda2_ad <<
" (AD) Error = " << std::setw(w)
128 << d2rda2 - d2rda2_ad << std::endl
129 <<
"d^2r/db^2 = " << std::setw(w) << d2rdb2 <<
" (analytic) == "
130 << std::setw(w) << d2rdb2_ad <<
" (AD) Error = " << std::setw(w)
131 << d2rdb2 - d2rdb2_ad << std::endl
132 <<
"d^2r/dadb = " << std::setw(w) << d2rdadb <<
" (analytic) == "
133 << std::setw(w) << d2rdadb_ad <<
" (AD) Error = " << std::setw(w)
134 << d2rdadb - d2rdadb_ad << std::endl
135 <<
"d^2r/dbda = " << std::setw(w) << d2rdadb <<
" (analytic) == "
136 << std::setw(w) << d2rdbda_ad <<
" (AD) Error = " << std::setw(w)
137 << d2rdadb - d2rdbda_ad << std::endl;
139 double tol = 1.0e-14;
140 if (std::fabs(r - r_ad) <
tol &&
141 std::fabs(drda - drda_ad) <
tol &&
142 std::fabs(drdb - drdb_ad) <
tol &&
143 std::fabs(d2rda2 - d2rda2_ad) <
tol &&
144 std::fabs(d2rdb2 - d2rdb2_ad) <
tol &&
145 std::fabs(d2rdadb - d2rdadb_ad) <
tol) {
146 std::cout <<
"\nExample passed!" << std::endl;
150 std::cout <<
"\nSomething is wrong, example failed!" << std::endl;
Sacado::Fad::DFad< double > DFadType
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
void func_deriv(double a, double b, double c, double &drda, double &drdb)
void func_deriv2(double a, double b, double c, double &d2rda2, double &d2rdb2, double &d2rdadb)
ScalarT func(const ScalarT &a, const ScalarT &b, const ScalarT &c)