lu_ratio.cpp¶

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;
}