Prev Next det_by_lu.hpp Headings

Source: det_by_lu
# ifndef CPPAD_DET_BY_LU_INCLUDED
# define CPPAD_DET_BY_LU_INCLUDED
# include <cppad/cppad.hpp>
# include <complex>

// BEGIN CppAD namespace
namespace CppAD {

// The AD complex case is used by examples by not used by speed tests 
// Must define a specializatgion of LeqZero,AbsGeq for the ADComplex case
typedef std::complex<double>     Complex;
typedef CppAD::AD<Complex>     ADComplex;
CPPAD_BOOL_UNARY(Complex,  LeqZero )
CPPAD_BOOL_BINARY(Complex, AbsGeq )

template <class Scalar>
class det_by_lu {
private:
     const size_t m;
     const size_t n;
     CppADvector<Scalar> A;
     CppADvector<Scalar> B;
     CppADvector<Scalar> X;
public:
     det_by_lu(size_t n_) : m(0), n(n_), A(n_ * n_)
     {    }

     template <class Vector>
     inline Scalar operator()(const Vector &x)
     {
          using CppAD::exp;

          Scalar         logdet;
          Scalar         det;
          int          signdet;
          size_t       i;

          // copy matrix so it is not overwritten
          for(i = 0; i < n * n; i++)
               A[i] = x[i];
 
          // comput log determinant
          signdet = CppAD::LuSolve(
               n, m, A, B, X, logdet);

# if 0
          // Do not do this for speed test because it makes floating 
          // point operation sequence very simple.
          if( signdet == 0 )
               det = 0;
          else det =  Scalar( signdet ) * exp( logdet );
# endif

          // convert to determinant
          det     = Scalar( signdet ) * exp( logdet ); 

# ifdef FADBAD
          // Fadbad requires tempories to be set to constants
          for(i = 0; i < n * n; i++)
               A[i] = 0;
# endif

          return det;
     }
};
} // END CppAD namespace
# endif

Input File: omh/det_by_lu_hpp.omh