Prev Next

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