det_by_minor.hpp

View page source

Source: det_by_minor

# ifndef CPPAD_DET_BY_MINOR_HPP
# define CPPAD_DET_BY_MINOR_HPP
# include <cppad/speed/det_of_minor.hpp>
# include <vector>

// BEGIN CppAD namespace
namespace CppAD {

template <class Scalar>
class det_by_minor {
private:
   //
   // m_
   // size for the matrix
   const size_t        m_;
   //
   // r_, c_
   // row and column indices so that minor is entire matrix.
   std::vector<size_t> r_;
   std::vector<size_t> c_;
   //
   // a_
   // temporary vector declared here to avoid reallocation for each use
   std::vector<Scalar> a_;
public:
   det_by_minor(size_t m) : m_(m) , r_(m + 1) , c_(m + 1), a_(m * m)
   {  //
      // r_, c_
      // values that correspond to entire matrix
      for(size_t i = 0; i < m; i++)
      {  r_[i] = i+1;
         c_[i] = i+1;
      }
      r_[m] = 0;
      c_[m] = 0;
   }
   //
   // operator()
   template <class Vector>
   Scalar operator()(const Vector &x)
   {  //
      // a_
      // copy from type Vector to std::vector<Scalar>
      for(size_t i = 0; i < m_ * m_; ++i)
         a_[i] = x[i];
      //
      // det
      // compute determinant of entire matrix
      Scalar det = det_of_minor(a_, m_, m_, r_, c_);
      //
# ifndef NDEBUG
      // r_, c_
      // values that correspond to entire matrix
      // (not const because det_of_minor uses r_, c_ for work space)
      for(size_t i = 0; i < m_; ++i)
      {  assert( r_[i] == i + 1 );
         assert( c_[i] == i + 1 );
      }
      assert( r_[m_] == 0 );
      assert( c_[m_] == 0 );
# endif
      return det;
   }

};

} // END CppAD namespace

# endif