|
|
Prev | Next |
int f2cad::dgefa_(doublereal *a, integer *lda, integer *n, integer *ipvt, integer *info);
int f2cad::dgesl_(doublereal *a, integer *lda, integer *n, integer *ipvt, doublereal *b, integer *job);
dgefa.f
and
dgesl.f
to solve the linear equation
\[
\left( \begin{array}{cc}
1 & 2 \\
3 & 4
\end{array} \right)
\left( \begin{array}{c}
f_0 \\
f_1
\end{array} \right)
=
\left( \begin{array}{c}
b_0 \\
b_1
\end{array} \right)
\]
Multiplying the top row by 3 and subtracting it from
the bottom row we obtain the equivalent equation:
\[
\left( \begin{array}{cc}
1 & 2 \\
0 & -2
\end{array} \right)
\left( \begin{array}{c}
f_0 \\
f_1
\end{array} \right)
=
\left( \begin{array}{c}
b_0 \\
b_1 - 3 b_0
\end{array} \right)
\]
It follows that
\[
\left( \begin{array}{c}
f_0 \\
f_1
\end{array} \right)
=
\left( \begin{array}{c}
b_1 - 2 b_0 \\
3 b_0 / 2 - b_1/ 2
\end{array} \right)
\]
The code below uses
the f2cad_link
routines
f2cad::Independent and f2cad::Dependent,
to define the function
\[
f(b) =
\left( \begin{array}{c}
b_1 - 2 b_0 \\
3 b_0 / 2 - b_1/ 2
\end{array} \right)
\]
We check that the derivative of this function,
calculated using the f2cad::Partial routine, satisfies
\[
f^{(1)} (b)
=
\left( \begin{array}{cc}
-2 & 1 \\
3/2 & -1 / 2
\end{array} \right)
\]
# include <f2cad/dgefa.hpp>
# include <f2cad/dgesl.hpp>
test_result dgesl(void)
{ bool ok = true;
// A is a 2 by 2 matrix in column major order
doublereal A[4];
A[0] = doublereal(1); A[2] = doublereal(2);
A[1] = doublereal(3); A[3] = doublereal(4);
// b is a column vector of length 2
doublereal b[2];
b[0] = doublereal(1);
b[1] = doublereal(2);
// declare independent variables
integer n = 2;
f2cad::Independent(n, b);
// f is a column vector of length 2
doublereal f[2];
integer i;
for(i = 0; i < n; i++)
f[i] = b[i];
// solve A * f = b
integer ipvt[2];
integer lda = n;
integer job = 0;
integer info;
f2cad::dgefa_ (A, &lda, &n, ipvt, &info);
f2cad::dgesl_ (A, &lda, &n, ipvt, f, &job);
// check info flag
ok &= (info == 0);
// construct the AD function object F : b -> f
// f[0] = b[1] - 2 b[0]
// f[1] = 3 b[0] / 2 - b[1]/ 2
integer m = 2;
f2cad::Dependent(m, f);
double p;
p = f2cad::Partial<doublereal>(0, 0);
ok &= f2cad::near_equal(p, -2., 1e-10, 1e-10);
p = f2cad::Partial<doublereal>(0, 1);
ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10);
p = f2cad::Partial<doublereal>(1, 0);
ok &= f2cad::near_equal(p, 3./2., 1e-10, 1e-10);
p = f2cad::Partial<doublereal>(1, 1);
ok &= f2cad::near_equal(p, -1./2., 1e-10, 1e-10);
if( ok )
return test_pass;
return test_fail;
}