Prev Next cppad_det_lu.cpp

CppAD Speed: Gradient of Determinant Using Lu Factorization

Operation Sequence
Note that the Lu factorization operation sequence depends on the matrix being factored. Hence we use a different ADFun object for each matrix.

compute_det_lu
Routine that computes the gradient of determinant using CppAD:
# include <cppad/vector.hpp>
# include <cppad/speed/det_by_lu.hpp>
# include <cppad/speed/uniform_01.hpp>

bool compute_det_lu(
     size_t                           size     , 
     size_t                           repeat   , 
     CppAD::vector<double>           &matrix   ,
     CppAD::vector<double>           &gradient )
{
     // -----------------------------------------------------
     // setup
     typedef CppAD::AD<double>           ADScalar; 
     typedef CppAD::vector<ADScalar>     ADVector; 
     CppAD::det_by_lu<ADScalar>          Det(size);

     size_t i;               // temporary index
     size_t m = 1;           // number of dependent variables
     size_t n = size * size; // number of independent variables
     ADVector   A(n);        // AD domain space vector
     ADVector   detA(m);     // AD range space vector
     
     // vectors of reverse mode weights 
     CppAD::vector<double> w(1);
     w[0] = 1.;

     // ------------------------------------------------------

     while(repeat--)
     {    // get the next matrix
          CppAD::uniform_01(n, matrix);
          for( i = 0; i < n; i++)
               A[i] = matrix[i];

          // declare independent variables
          Independent(A);

          // AD computation of the determinant
          detA[0] = Det(A);

          // create function object f : A -> detA
          CppAD::ADFun<double> f(A, detA);

          // evaluate and return gradient using reverse mode
          gradient = f.Reverse(1, w);
     }
     return true;
}

Input File: speed/cppad/det_lu.cpp