Prev Next Index-> contents reference index search external Up-> CppAD Appendix numeric_ad LuRatio lu_ratio.cpp Appendix-> Faq directory Theory glossary Bib wish_list whats_new deprecated compare_c numeric_ad addon License numeric_ad-> BenderQuad opt_val_hes LuRatio LuRatio-> lu_ratio.cpp lu_ratio.cpp Headings

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
LuRatio: Example and Test
# include <cstdlib> // for rand function # include <cassert> # include <cppad/cppad.hpp> namespace { // Begin empty namespace CppAD::ADFun<double> *NewFactor( size_t n , const CPPAD_TESTVECTOR(double) &x , bool &ok , CPPAD_TESTVECTOR(size_t) &ip , CPPAD_TESTVECTOR(size_t) &jp ) { using CppAD::AD; using CppAD::ADFun; size_t i, j, k; // values for independent and dependent variables CPPAD_TESTVECTOR(AD<double>) Y(n*n+1), X(n*n); // values for the LU factor CPPAD_TESTVECTOR(AD<double>) LU(n*n); // record the LU factorization corresponding to this value of x AD<double> Ratio; for(k = 0; k < n*n; k++) X[k] = x[k]; Independent(X); for(k = 0; k < n*n; k++) LU[k] = X[k]; CppAD::LuRatio(ip, jp, LU, Ratio); for(k = 0; k < n*n; k++) Y[k] = LU[k]; Y[n*n] = Ratio; // use a function pointer so can return ADFun object ADFun<double> *FunPtr = new ADFun<double>(X, Y); // check value of ratio during recording ok &= (Ratio == 1.); // check that ip and jp are permutations of the indices 0, ... , n-1 for(i = 0; i < n; i++) { ok &= (ip[i] < n); ok &= (jp[i] < n); for(j = 0; j < n; j++) { if( i != j ) { ok &= (ip[i] != ip[j]); ok &= (jp[i] != jp[j]); } } } return FunPtr; } bool CheckLuFactor( size_t n , const CPPAD_TESTVECTOR(double) &x , const CPPAD_TESTVECTOR(double) &y , const CPPAD_TESTVECTOR(size_t) &ip , const CPPAD_TESTVECTOR(size_t) &jp ) { bool ok = true; using CppAD::NearEqual; double eps99 = 99.0 * std::numeric_limits<double>::epsilon(); double sum; // element of L * U double pij; // element of permuted x size_t i, j, k; // temporary indices // L and U factors CPPAD_TESTVECTOR(double) L(n*n), U(n*n); // Extract L from LU factorization for(i = 0; i < n; i++) { // elements along and below the diagonal for(j = 0; j <= i; j++) L[i * n + j] = y[ ip[i] * n + jp[j] ]; // elements above the diagonal for(j = i+1; j < n; j++) L[i * n + j] = 0.; } // Extract U from LU factorization for(i = 0; i < n; i++) { // elements below the diagonal for(j = 0; j < i; j++) U[i * n + j] = 0.; // elements along the diagonal U[i * n + i] = 1.; // elements above the diagonal for(j = i+1; j < n; j++) U[i * n + j] = y[ ip[i] * n + jp[j] ]; } // Compute L * U for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { // compute element (i,j) entry in L * U sum = 0.; for(k = 0; k < n; k++) sum += L[i * n + k] * U[k * n + j]; // element (i,j) in permuted version of A pij = x[ ip[i] * n + jp[j] ]; // compare ok &= NearEqual(pij, sum, eps99, eps99); } } return ok; } } // end Empty namespace bool LuRatio(void) { bool ok = true; size_t n = 2; // number rows in A double ratio; // values for independent and dependent variables CPPAD_TESTVECTOR(double) x(n*n), y(n*n+1); // pivot vectors CPPAD_TESTVECTOR(size_t) ip(n), jp(n); // set x equal to the identity matrix x[0] = 1.; x[1] = 0; x[2] = 0.; x[3] = 1.; // create a fnction object corresponding to this value of x CppAD::ADFun<double> *FunPtr = NewFactor(n, x, ok, ip, jp); // use function object to factor matrix y = FunPtr->Forward(0, x); ratio = y[n*n]; ok &= (ratio == 1.); ok &= CheckLuFactor(n, x, y, ip, jp); // set x so that the pivot ratio will be infinite x[0] = 0. ; x[1] = 1.; x[2] = 1. ; x[3] = 0.; // try to use old function pointer to factor matrix y = FunPtr->Forward(0, x); ratio = y[n*n]; // check to see if we need to refactor matrix ok &= (ratio > 10.); if( ratio > 10. ) { delete FunPtr; // to avoid a memory leak FunPtr = NewFactor(n, x, ok, ip, jp); } // now we can use the function object to factor matrix y = FunPtr->Forward(0, x); ratio = y[n*n]; ok &= (ratio == 1.); ok &= CheckLuFactor(n, x, y, ip, jp); delete FunPtr; // avoid memory leak return ok; } 
Input File: example/general/lu_ratio.cpp