/* ---------------------------------------------------------------------------- f2cad: Convert Fortran to C++ AD Types: Copyright (C) 2005-10 Bradley M. Bell Authors: Bradley M. Bell bradbell at washington dot edu License: GNU Public License Version 2 or higher -----------------------------------------------------------------------------*/ /* $begin dgesl$$spell hpp blas df Linpack dgefa dgesl bool doublereal ipvt lda$$$section dgesl$$index blas, dgesl$$ $index dgesl, blas$$head Prototype$$$include p.omh/dgefa.p.omh$$include p.omh/dgesl.p.omh$$ $children% omh/dgesl.f.omh %$$head Fortran Source$$$cref dgesl.f$$head Description$$ This example uses the Linpack routines $codei%dgefa.f%$$and codei%dgesl.f%$$ to solve the linear equation$latex $\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: latex $\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 $latex $\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 xref/f2cad_link/$$ routines$code f2cad::Independent$$and code f2cad::Dependent$$, to define the function $latex $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 code f2cad::Partial$$ routine, satisfies$latex $f^{(1)} (b) = \left( \begin{array}{cc} -2 & 1 \\ 3/2 & -1 / 2 \end{array} \right)$ $$codep */ # include # include bool 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(0, 0); ok &= f2cad::near_equal(p, -2., 1e-10, 1e-10); p = f2cad::Partial(0, 1); ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10); p = f2cad::Partial(1, 0); ok &= f2cad::near_equal(p, 3./2., 1e-10, 1e-10); p = f2cad::Partial(1, 1); ok &= f2cad::near_equal(p, -1./2., 1e-10, 1e-10); return ok; } /*$$ \$end */