Prev Next Index-> contents reference index search external Up-> f2cad library dgefa f2cad-> license Install get_started prototype library run.cpp utility add2lib.sh whats_new library-> d9knus d9lgmc daxpy dbesks dbskes dcsevl ddot dgamlm dgamma dgefa dgesl dint dscal idamax initds xerror dgefa-> dgefa.f Headings-> Prototype Fortran Source Description

dgefa

Prototype
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);

Fortran Source
dgefa.f

Description
This example uses the Linpack routines dgefa.f and dgesl.f to solve the linear equation   $\left( \begin{array}{cc} 1 & 2 \\ 0 & 3 \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)$  Using the f2cad_link routines f2cad::Independent and f2cad::Dependent, this defines the function   $f(b) = \left( \begin{array}{c} b_0 - 2 b_1 / 3 \\ b_1 / 3 \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} 1 & - 2 / 3 \\ 0 & 1 / 3 \end{array} \right)$   # include <f2cad/dgefa.hpp> # include <f2cad/dgesl.hpp> test_result dgefa(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(0); A[3] = doublereal(3); // 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[1] = b[1] / 3 // f[0] = b[0] - 2 * b[1] / 3 integer m = 2; f2cad::Dependent(m, f); double p; p = f2cad::Partial<doublereal>(0, 0); ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 1); ok &= f2cad::near_equal(p, -2./3., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 0); ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 1); ok &= f2cad::near_equal(p, 1./3., 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dgefa.cpp