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